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Chapter  1 
INTRODUCTION 

The  problems  of  linear  system  synthesis  and  identification  are 
closely  related.  In  system  synthesis,  one  tries  to  synthesise  an 
input-output  relation  that  will  yield  a specified  behavior.  The 
first  step  usually  consists  of  determining  either  the  impulse  re- 
sponse or  transfer  function  of  the  system.  The  impulse  response  is 
then  approximated  by  a sum  of  exponentials  while  the  transfer  func- 
tion is  approximated  by  a realizable  rational  function.  Since  the 
physical  system  does  not  exist  during  the  synthesis  phase,  system 
synthesis  is  characterized  by  the  feature  that  one  deals  with  ideal- 
ized analytical  functions. 

This  is  not  the  case  with  system  identification.  Here  the  system 
actually  exists  and  the  identification  process  is  based  upon  measured 
data  at  the  system  input  and  output  ports.  This  data  is  usually  con- 
taminated by  measurement  errors,  thermal  noise,  and  spurious  signals. 

t 

In  addition,  the  record  lengths  are,  by  necessity,  finite. 

In  system  identification,  as  in  system  synthesis,  complex  expo- 
nential functions  provide  a natural  basis  for  the  impulse  response. 

However,  in  some  cases  the  impulse  response  cannot  be  exactly  repre- 
sented as  a sum  of  exponentials.  The  identification  problem  then 


becomes  an  approximation  problem. 
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The  intrinsic  difference  between  synthesis  and  identification  is 
that  the  former  seeks  a system,  not  yet  in  existence,  to  approximate 
an  operator  which  relates  idealized  inputs  and  outputs,  while  the  latter 
attempts  to  describe  an  existing  unknown  system  through  measured  input- 
output  records.  Since  both  problems  emphasize  the  input-output  behavior 
of  an  unknown  system,  the  strategies  involved  with  their  solution  are 
similar.  However,  they  are  distinguished  on  the  basis  of  the  type  of 
data  provided.  Whereas  a synthesis  problem  is  characterized  by  syn- 
thetic data,  an  identification  problem  is  characterized  by  measured 
data  records. 

The  goal  of  this  dissertation  is  to  develop  a general  parametric 
approximation/ identification  technique  for  linear  time-invariant  systems. 
The  technique  models  the  system  impulse  response  as  a sum  of  exponentials. 
It  is  shown  that  this  technique  can  be  extended  to  slowly  time-variant 
and  weakly  nonlinear  systems. 

Chapter  1 presents  some  useful  definitions  and  discusses  system 
classification  with  respect  to  input-output  relations.  Beginning  with 
the  general  class  of  systems,  attention  is  first  focussed  on  linear 
time-invariant  systems  and  then  special  ized  further  to  the  subset  of 
rational  systems  within  this  class.  The  suitability  of  the  exponential 
family  of  functions  for  modeling  an  impulse  response  is  examined.  The 


1.1  System  Classification  and  Definitions 

Before  addressing  the  system  identification  problem,  a brief 
discussion  about  the  classification  of  systems,  their  pertinent  quali- 
tative properties,  and  some  definitions  are  presented. 

In  general,  characterizat  on  of  the  input-output  properties  of 
a system  implies  an  operator  H that  operates  on  a class  of  signals 
(functions)  {x}  and  results  in  another  class  of  signals  {y}.  A system 
is  relaxed  (or  at  rest)  at  time  t = t if  inputs  are  not  applied  over 

-°°  < t < t and/or  tne  output  signal  depends  solely  upon  the  operation 
— — o 

of  H on  the  input  signal  over  tQ  < t < 00 . In  the  most  general  case, 
an  operator  H might  involve  memory  and  be  nonlinear,  time-variant, 
anticipative , and  ambiguous. 

The  above  properties  are  defined  as  follows  [1]: 


1-4 


for  any  input  signals  x^,  and  for  ar.y  real  numbers  i = l,2,...,n. 

Otherwise  the  relaxed  system  is  nonlinear. 

(c)  Let  Q denote  the  time-shift  operator  with  shift  a.  A relaxed  system 
a 

is  time  invariant  (stationary  or  fixed)  if  and  only  if 


HQ  x = Q Hx 

a a 


(1-1-2) 


for  any  input  signal  x,  and  for  any  real  number  a.  Otherwise  the 
system  is  time-variant. 

(d)  A system  is  causal  (nonanticipatory)  if  and  only  if  the  output 
signal  y(t)  depends  solely  upon  the  operation  of  H on  the  present 
and  past  of  the  incoming  signal  x(t).  If  y(t)  also  depends  upon 
future  values  of  x(t),  the  system  is  anticipatory. 

(e)  A system  is  said  to  be  unique  if  and  only  if  there  exists  a single 
output  signal  y(t)  for  any  input  signal  x(t).  Otherwise,  the  system 
is  ambiguous. 

In  this  dissertation  only  physical  systems  which  are  causal,  unique, 
slowlv  time  - variant»  and  weakly  nonlinear  with  memory  are  considered. 

In  general,  the  input-output  relationship  for  a weakly  nonlinear  system 
can  be  characterized  in  terms  of  the  Volterra  representation  which  is 
expressed  as  [2] 


y(t)  - l ...  h (t  I ) n x(t,T  )dl 

i=l  J 1 1=1  3 3 * 

— OO  —00  J 


(1-1-3) 


This  represents. ion  was  first  utilized  by  Norbert  Wiener  [3]  in  the  analy- 
sis of  nonlinear  systems.  For  physical  systems  which  are  initially 
relaxed  at  time  tQ,  the  input-output  relation  can  be  written  as 
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y(t) 


i 


t-t 
r o 


(t^ , • 


i 

,T  ) n x(t  ,T  )dx  , t > t 
1 j=l  * i ~ ° 


(1-1-4) 


In  addition,  the  infinite  sum  in  (1-1-4)  can  be  adequately  approximated 
by  a finite  number  of  terms  for  systems  which  are  weakly  nonlinear.  Also, 
for  systems  which  are  time-invariant,  (1-1-4)  becomes 


t-t  t-t 
o o 


y(t)  = I 

i=l 


o o 

~T 


hi(T1,...,Ti)  n x(t  - T.)dT  , t _>  t (1-1 
1=1  J J 


-1-5) 


Alternatively,  (1-1-5)  can  be  written  as 
t t 


N 


y(t)  = l 

i=l  J 


t t 
. o o 


h (t-T.  , . . . ,t-T.)  IT  x(T.)dT.,  t > t 

11  i i=1  1 j - O 


(1-1-6) 


In  practice,  there  exists  a broad  class  of  electronic  circuits 
which  are  weakly  nonlinear.  Values  of  N from  3 to  5 are  commonly  used 
when  modeling  amplifiers  in  communications  receivers  [4],  For  N=l, 
(1-1-6)  describes  the  input-output  relationship  of  a linear,  time- 
invariant,  causal  system  with  memory  which  is  relaxed  at  time  tQ. 


1.2  Input-output  Stability  for  Linear  Systems 

One  of  the  important  characteristics  of  a physical  system  is  its 
stability.  Since  an  identification  process  utilizes  information  pro- 
vided by  the  input  and  output  signals  of  the  system,  it  is  convenient 
to  describe  a system's  stability  in  terms  of  its  impulse  response. 


J 
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The  following  useful  theorems  [1]  are  referred  to  in  the  later 
development : 

Theorem  1-2-1 

A relaxed  system  at  t = t is  bounded-input  bounded  output  stable 
(BIBO)  if  and  only  if  there  exists  a finite  number  K such  that 

t 

J |h(t,T) |dx  £ K < ~ (1-2-1) 

t 

o 

for  all  t < t < °°. 
o — — 

For  the  time- invariant  case,  the  theorem  becomes  stronger. 


Theorem  1-2-2 

Consider  a linear  time- invariant  system  relaxed  at  t • t , which 
is  described  by  the  following  relation 


If 


t 


y(t)  = 


J 

t 

o 


h (t— t)x (t) dt  . 


oo 


|h(t)  | dt  £ K 


< 


CO 


(1-2-2) 


(1-2-3) 


for  some  constant  K,  then  we  have  the  following: 

(a)  If  x is  a periodic  function  with  period  T for  all  t _>  tQ, 
then  the  output  y tends  to  a periodic  function  with  the  same 
period  (not  necessarily  of  the  same  waveform). 

(b)  If  x is  bounded  and  tends  to  a constant,  then  the  output  y is 


bounded  and  tends  to  a constant. 
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(c)  If  x is  of  finite  energy,  that  is, 

.00 

(j  |x(t)|2dt)1/2  < Kx  < » , 


(1-2-4) 


then  the  output  is  also  of  finite  energy;  that  is,  there  exists  a finite 
that  depends  on  K,  such  that 


(j  |y(t) j2dt)1/2  < K2  < °° 


(1-2-5) 


(d)  The  system  is  BIBO  stable  (This  is  a corollary  of  Theorem  1-2-1). 

For  systems  described  by  proper  rational  transfer  functions  the  following 
theorem  applies. 

Theorem  1-2-3 

A system  relaxed  at  t = tQ  that  is  described  by  a proper  rational 
function  H(S)  (or  equivalently,  an  impulse  response  h(t)  which  is  a 
sum  of  complex  exponentials)  Is  BIBO  stable  if  and  only  if  all  the  poles 
of  H(S)  are  in  the  open  left  half  S-plane  (i.e.,  excluding  the  imaginary 
axis).  All  the  poles  of  H(S)  then  have  negative  real  parts. 

The  only  requirement  imposed  on  the  kernel  h(')  in  theorems  1-2-1 
and  1-2-2  is  that  it  should  be  absolutely  integrable.  This  fact  does 
not  imply  boundedness  or  that  h(t)  approaches  zero  as  t ■*  °°.  On  the 
other  hand,  the  impulse  response  of  theorem  1-2-3  is  bounded  and  ap- 
proaches zero  as  t + ",  We  may  write  for  this  case 


k mi 


MO  -II  A.  .t*-1.  1W,  t > 0 
< i-1  j-1  1J 


; 0 > t 


(1-2-6) 
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where  A + 0,  a.  - Real  [s.]  < 0,  s.  1 s.  for  1,  ^ i„.  m.  denotes 
lmi  1 1 i1  12  1 2 i 

th«*  order  of  the  pole  at  s.,  The  total  number  of  poles  in  h(t)  is 
k 1 

. iven  by  J m.  = n.  Using  (1-2-6)  we  write 
1-1  1 


k mi 


i m . 

. , st  k i at 

l l ! I |A„U3-le  1 


h(t) 


i-1  j-1 


ij 


i-1  j-1 


ij 


k t 

£ £ max 

i=l 


'il1 


d-2-7) 


Inequality  (1-2-7)  shows  the  boundedness  of  h(t).  It  is  clear  that 

lim  h(t)  * 0 since  0^  < 0. 

t-*00 

As  a matter  of  fact,  (1-2-6)  not  only  belongs  to  L^tO,00]  space 

p 

(i.e.tis  absolutely  integrable  over  [0,°°]),  but  also  belongs  to  L [O,00] 
where  P _>  1,  since 


r k i , . 8 t.‘ 

I l l A tJ"Ae  1 | dt  < co  . (1-2-8) 

b 1=1  j=1 


In  L [0,°°]  (Hilbert  space)  the  infinite  sum 
„ m , 


I I 

i-1  j-1 


a J-1  8it 

V 6 


(1-2-9) 


provides  a complete  representation  for  piecewise  continuous  square 
integrable  functions,  assuming  proper  choices  of  the  poles  {s^}.  As  a 
matter  of  fact,  many  choices  of  {s  } exist.  From  Szasz’s  form  of 

i-i  V 

Muntz's  theorem  [5],  all  that  is  required  for  the  set  {tJ  e } to 

2 

be  complete  in  L [0,°°]  is  that 
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00  Real  [s.] 

1  

1=1  1 + |si  + \\2 


(1-2-10) 


For  example,  possible  complete  sets  are 
(a)  {e  t,  e 2t,...,e  it:  ...} 


r -(l+^T)t  -2(l+vCT)t  -i(l+»£T)t 

*- e >6  9 • * 9 e f • ■ 


m1=  j = i 

i = 1,2,.. 


In  practice,  when  approximating  a waveform  by  (1-2-9),  a finite  set 
is  used  and  an  error  is  involved  in  the  representation  which  one  tries 
to  minimize  by  a proper  choice  of  {s^}..  Hence,  representations  in 
the  form  of  (1-2-6)  result. 

It  can  be  shown  [6]  that  a function  which  both  belongs  to 

1 9 

L [0,°°]  and  is  bounded  also  belongs  to  L‘  [0,°°].  In  other  words,  the 

boundedness  constraint  added  to  an  function  is  sufficient  to 

2 

guarantee  membership  in  L . (In  fact,  membership  is  also  guaranteed 
in  LP,  P _>  1). 

In  physical  systems,  the  boundedness  requirement  is  not  a severe 
constraint  (i.e.,  physical  stable  systems  are  inherently  bounded.; 

It  is  true  that  theorem  1-2-1  does  allow  h(*)  to  contain  a countab? j 
linear  combination  of  impulses  provided  their  coefficients  constitute 
an  absolutely  converging  series.  However,  as  is  well  known,  there 
is  no  real  system  for  which  the  impulse  response  is  an  ideal  impulse. 
Therefore,  all  transfer  functions  belonging  to  physical  systems  have 
a zero  at  °°.  However,  transfer  functions  cannot  be  identically  zero 


over  any  nonzero  frequency  interval.  The  Paley-Wiener  criterion  [7] 
allows  at  most  a countable  number  of  zeros  to  exist  on  the  imaginary 
axis  in  the  s plane  for  any  transfer  function  related  to  a causal 
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L 


impulse  response. 


1.3  The  Rational  Linear  System 

A general  rational  linear  system  is  described  in  the  time  domain 
by  its  impulse  response  as  given  in  (1-2-6).  The  Laplace  transform  of 
(1-2-6)  yields  the  rational  transfer  function 

k "i  lv1 


H(S)  - l l *,, 

J 1 1 1 


« ^ 4-1  A n 

i“l  J-l  - (S  " J d SJ 

1=0  j 


(1-3-1) 


where  m < n.  This  is  a frequency-domain  description.  Clearly,  (1-3-1) 
has  a countable  number  of  finite  zeros  (exactly  m)  and  a zero  at 
infinity. 

One  way  to  describe  the  input-output  relations  for  such  a system 
is  through  the  convolution  integral.  Using  (1-2-6)  in  (1-2-2),  there 
results 


k mi 


y(t)  -l  l A 


i=l  j-1 


ij 


Vt"0 

(t-T)-1  "e  x(t)dl 


(1-3-2) 


where  the  system  is  assumed  to  be  relaxed  at  t = 0.  Taking  the  Laplace 
transform  of  (1-3-2),  we  obtain 

m . 

I V 

I A 


Y (S)  - H(S)X(S)  - 


i-0 


n . 

I DjsJ 

j-0  3 


X(S) 


(1-3-3) 


Premultiplying  (1-3-3)  by  the  denominator  polynomial  yields 


l D.SJY(S)  - l N1S1X(S)  . 

j-0  J 


(1-3-4) 


i-0 
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The  inverse  Laplace  transform  of  (1-3-4)  is 


r (i)  v ( 

T D y ' (t)  = l N.  x 

j =0  J i=0  1 


(1-3-5) 


where  the  ini  ial  conditions  are  given  by  yv  ' (0)  = 0,  k = 0,1,2, ... ,n-l 
x(e)(0'  = 0,  i = 0,1,2,. ..,m-l. 

Observe  that  y^-^(t)  is  the  jth  derivative  of  y(t)  and  x^^(t) 
is  the  ith  derivative  of  x(t).  Thus,  (1-3-5)  is  a differential  equa- 
tion which  represents  the  input-output  relation  of  a causal,  time- 

invariant,  linear  system.  For  the  more  general  situation  in  which  the 

00 

initial  conditions  are  nonzero  (i.e.,  y (0)  y*  0,  k * 0,1,..., n-1  and 

(£) 

x (0)  y*  0,  Z = 0,1,2, ... ,m-l) , the  Laplace  transform  of  (1-3-5)  is 

m . n-1 

l Nisl  l vq 

Y(S)  = X(S)  ^ ^ , (1-3-6) 

n . n . 


I D SJ  £ DiSJ 

j=0  J j=0  J 


where 


~l  C Sq  = l J[1  D Sky(j"k_1)(0)  - l l N S£x(i"<l"1)(0), 

q=0  q j=l  k=0  J i=l  l-  0 


l D y(J-q  ] )(0)  - l N1xU-q-i;  (0). 


q L 1 

q j=q+i  J 


i*  q+1 


If  (1-3-1)  is  irreducible,  then  (1-3-6)  is  said  to  represent  the  Laplace 
transform  of  the  differential  equation  of  an  observable  and  controllable 
system  [1],[8].  Only  irreducible  transfer  functions  are  considered  in 
this  dissertation. 
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1.4  Literature  Review  and  Problem  Statement 

Many  investigators  have  been  Involved  with  various  aspects  of 
the  system  identification  problem  since  the  late  1940' s.  The  proposed 
techniques  can  be  partitioned  between  those  that  utilize  parametric 
versus  nonparametric  models.  Nonparametric  models  are  characterized 
by  the  fact  that  the  linear  impulse  response  and/or  transfer  function 
is  presented  in  either  tabular  or  graphical  form.  On  the  other  hand, 
parametric  models  result  in  a numerical  set  of  values  for  parameters 
associated  with  a set  of  basis  functions  used  to  model  the  impulse 
response  and/or  transfer  function.  This  dissertation  is  devoted  to 
the  parametric  approach  in  which  the  basis  functions  for  the  impulse 
response  are  chosen  to  be  a set  of  complex  exponential  time  functions 
where  the  coefficients  and  exponents  are  the  parameters  to  be  deter- 
mined. 


Historically,  the  first  attempts  at  the  identification  problem 
were  linked  to  system  synthesis  where  exponential  approximations  were 
generated  to  represent  a prescribed  impulse  response.  The  simplest 
solutions  were  obtained  by  first  specifying  the  exponents  a priori  and 
then  solving  a linear  set  of  equations  for  the  coefficients  based  upon 
either  linear  interpolation,  the  least-square  approach,  or  minimization 
of  the  integrated  squared  error  [9],  By  specifying  the  exponents  a 
priori,  the  difficult  part  of  the  problem  was  avoided.  However,  even 
when  the  exponents  are  not  specified  a priori,  the  coefficients  may 
be  determined  as  mentioned  above.  Appendix  A provides  a tutorial  review 
for  the  rase  in  whirh  the  set  of  hasis  functions  is  chosen  in  advance. 
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Other  approximations  for  the  impulse  response  involved  known  ortho- 
normal sets  of  exponential  basis  functions  [10],  [11],  [12], 

Since  an  a priori  choice  of  exponents  is  uneconomical  (i.e.,  the 
number  of  exponents  needed  for  an  adequate  approximation  may  be  un- 
necessarily large),  attention  was  focused  on  the  nonlinear  problem  of 
obtaining  the  "best"  exponents  according  to  some  criterion.  One  pro- 
cedure was  suggested  by  Prony  [13]  as  early  as  1795.  This  method  uses 
a set  of  2N  equidistant  data  points  to  generate  2N  nonlinear  equa- 
tions in  the  N unknown  exponents  and  their  corresponding  N unknown 
coefficients.  In  an  effort  to  obtain  increased  accuracy  and  noise 
immunity,  Prony's  method  was  extended  to  make  use  of  more  than  2N 
data  points  in  the  search  for  the  N exponentials  [14].  Other  forms 
of  the  overdetermined  Prory  method  were  used  in  the  areas  of  speech 
analysis  [15]  and  electromagnetic  fields  [16],  [17],  [18].  All  of  the 
above  techniques  result  in  an  estimate  of  the  system  characteristic 
equation  which  is  solved  for  the  unknown  exponents. 

Other  methods,  arising  from  the  system  synthesis  problem,  involve 
iterative  schemes  to  determine  the  exponents  and  residues  [19],  [20], 
[21],  [22].  These  methods  attempt  to  fit  a closed  form  expression 
for  the  impulse  response  by  a sum  of  exponentials.  The  procedure  is 
optimized  by  minimizing  either  the  integrated-squared  error  or  least- 
squares  error. 

Whereas  Prony's  method  uses  delay  operators  to  obtain  the  charac- 


teristic equation  of  a difference  equation,  other  techniques  have 
used  either  differential  or  integral  operators  to  obtain  the  charac- 
teristic equation  of  a differential  equation  [23],  [24],  [25],  [26], 
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As  before,  the  exponents  are  determined  from  the  characteristic  equa- 
tion. 

All  of  the  previous  techniques  can  be  used  in  the  identification 
problem  only  when  the  waveform  being  approximated  is  the  impulse  response 
of  the  system.  This  implies  that  the  input  is  restricted  to  be  an 
impulse.  However,  other  investigators  have  shown  that  the  identifica- 
tion problem  can  be  solved,  given  the  input  and  output  signals,  even 
when  the  input  has  some  arbitrary  waveform.  Gradient  iterative  tech- 
niques using  general  inputs  are  discussed  in  [27]  and  [28],  An  exten- 
sion of  Prony's  method  for  arbitrary  inputs  is  presented  in  [8]. 

Finally,  the  approaches  given  in  [25]  and  [26]  are  generalized  for 
general  input  signals  in  [29],  [30],  and  [31]. 

As  is  evident  from  the  survey,  many  different  approaches  have  been 
proposed  for  the  linear  system  identification  problem.  Each  approach 
has  certain  advantages  and  disadvantages.  In  addition,  the  various 
techniques  have  been  proposed  without  insight  into  how  they  are  related. 
In  this  dissertation,  a generalized  parametric  model  is  presented  which 
unifies  those  techniques  that  generate  the  system  characteristic  equa- 
tion. Each  of  these  techniques  results  by  choosing  the  appropriate 
operator  in  the  general  model.  These  techniques  all  rely  upon  the 
inversion  of  a Hilbert  type  of  inner-product  matrix  which  is  basically 
ill-posed.  One  advantage  of  the  generalized  model  is  that  it  enables 
the  flexibility  in  a specific  problem  to  select  an  operator  which  tends 
to  diagonalize  the  matrix.  This  is  of  considerable  help  when  inverting 
the  matrix.  The  generalized  model  also  enables  the  choice  of  operators 
for  the  purpose  of  increasing  noise  immunity.  Some  of  the  techniques 


1-15 


proposed  in  the  literature  required  the  availability  of  data  over  a 
semi-infinite  interval.  In  practical  applications  the  data  is  truncated. 
This  results  in  a truncation  error.  One  of  the  by-products  of  the 
general  model  is  the  ability  to  deal  with  portions  of  a data  record 
without  incurring  truncation  errors.  In  addition,  most  of  the  itera- 
tive schemes  proposed  in  the  literature  suffer  from  either  convergence 
problems  or  the  lack  of  a proof  of  convergence.  The  generalized  model 
may  be  used  in  conjunction  with  an  iterative  scheme  by  providing  a 
good  first  guess  at  the  solution.  This  eases  the  convergence  problem. 
Finally,  the  generalized  model  is  readily  extended  to  handle  complex 
data  as  opposed  to  real  data. 

1 . 5 Dissertation  Outline 

Chapter  2 is  devoted  to  the  development  of  the  gei:  ral  parametric 
model  for  identification.  Since  the  exponential  family  of  functions 
plays  a major  role  in  the  suggested  model,  we  examine  in  Section  2.1 
the  behavior  of  this  family  when  subjected  to  linear  operators.  Also, 
in  2.2,  their  correlation  properties  are  investigated.  In  2.3  the 
general  parametric  model  is  actually  developed,  and  in  2.4  the  identi- 
fication procedure  is  demonstrated  through  examples  where  a variety  of 
operators  are  used  and  their  properties  examined. 

Chapter  3 introduces  some  supplemental  tools.  First  the  concept 
of  the  S-plane  "growing  tree"  is  introduced,  its  properties  and  structure 
discussed,  and  potential  benefits  when  applied  to  synthesis  and  identi- 
fication problems  pointed  out. 


The  need  for  a good  initial  guess  in  iterative  techniques  makes 
the  model  discussed  in  chapter  2 an  a' tractive  candidate  as  a generator 
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of  such  a guess.  The  use  of  such  a scheme  to  obtain  the  initial  guess 
is  demonstrated  through  an  example. 

Chapter  4 treats  the  identification  of  a weakly  nonlinear  system. 
The  system  is  modeled  through  a Volterra  series  representation.  An 
identification  procedure  for  the  higher  order  impulse  responses  is 
outlined  and  demonstrated  through  an  example. 

A summary  of  the  significant  results  and  conclusions  are  presented 
in  chapter  5. 

Appendix  A presents  some  tutorial  material  about  possible  approxi- 
mation schemes. 
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Chapter  2 

PARAMETRIC  IDENTIFICATION  SCHEMES 
VIA  LINEAR  OPERATORS 

The  general  parametric  identification  model  is  developed  in  this 
chapter.  In  this  model  it  is  assumed  that  the  data  available  consists 
of  input  and  output  records  solely  (i.e.,  black  box  identification).  The 
model  requires  that  those  records  be  processed  by  linear  operators.  The 
exponential  family  of  functions,  being  a natural  candidate,  is  the  one 
chosen  to  model  the  identified  system.  The  behavior  of  that  family,  when 
subjected  to  linear  operations,  is  first  investigated.  Also,  the  correla- 
tion properties  of  this  family  of  functions  is  examined. 

It  is  shown  that  the  group  of  operators  that  can  be  used  for  identi- 
fication purposes  is  large.  The  ease  of  use  of  a certain  operator  depends 
upon  the  ease  of  obtaining  the  parameter  transformation  associated  with 
the  particular  operator.  Different  operators  will  affect  differently  a 
given  identification  situation.  In  general,  one  should  choose  an  operator 
that  is  "tailored"  for  the  identification  problem  at  hand.  Some  operators 
are  inherently  ill-posed  when  used  in  any  identification  scheme  outlined 
by  the  general  model. 

Several  examples  are  presented  using  a variety  of  operators  and 
different  approximation  schemes  from  Appendix  A. 


2 . 1 Effect  of  Linear  Operator  on  the  Exponential  Basis 

t h 

In  general,  an  n -order  complex  exponential  basis,  with  negative 

real  parts  for  the  complex  frequencies,  is  presented  by 

i = 1,2, . . . ,k  k 

„ . . , „ y m.  = u 

*1  S t j = 1,2,. ..  ,m  t l 

(tJ  ‘ >>  real  [,,1  < 0*  (2-1‘1> 

t > 0. 


A linear  combination  representing  either  a function  which  belongs 
to  the  space  spanned  by  this  basis  or  the  projection  onto  the  space  of 
a function  which  does  not  belong  to  the  space  is  given  by 


k m.  St 

e 


g(o  = i r A.jt 


i=l  j=l 


t > 0 


(2-1-2) 


In  this  section  the  behavior  of  such  representations  when  subjected  to 
certain  operators  is  investigated.  The  first  operator  to  be  considered 
is  the  differential  operator.  Operating  upon  (2-1-2)  with  a differen- 
tiator (for  t^O),  there  results 


where 


g 


at  i=l  j*l  1 J 


- Vij  + jAi(j+i) 

^lim,  = Vim, 


J“l,2,. . . , (m^-1) 


(2-1-3) 


Subjecting  (2-1-3)  to  a second  differentiator,  we  obtain 

t > 0 


m dg  (t)  k m S.t 

'"’w-i-,  i • 1 

1=1  j=l  J 


(2-1-4) 


where 


A- . . = S .A  . . + jA  . , . _ , 
2ij  l lij  J H(j+1) 


A0 . — S .A  . 

2im.  l Inn. 

l l 


j=l,2,...,(m.-l) 


In  general,  after  (2-1-2)  has  been  subjected  to  £ differentiations, 

g()(t)=--S-— = l l A tJ_  e 1 ; t>0  (2-1-5) 

d i=l  j=l 


where  A...  = S.A..  + jA,„  . 

£ij  i (£-l)ij  (£-1) i (j+1) 


£im.  ^i^(£-l)im. 
i l 


j=l,2, . . . , (m^-1) . 


Looking  at  (2-1-5),  we  observe: 

(a)  Equation  (2-1-2)  is  infinitely  differentiable  for  all  t ^ 0. 

(£) 

(b)  g (t)  is  a linear  combination  using  the  original  basis  (2-1-1) 
for  any  £ (i.e.,  space  invariance  is  experienced  under  repeated 
operation  of  the  differential  operator  on  functions  in  the  space 
spanned  by  (2-1-1)). 

(c)  For  any  integer  value  of  £^,  the  basis  set 


k m. 

I l I1  A 
i-1  j-1 


Uj 


i-1  ^i*" 

tJ  e };  £ = £1,(£1+1) (^4-n-l)  (2-1-6) 


spans  exactly  the  same  space  as  does  the  set  (2-1-1). 

(d)  Any  m+1  set  (i.e.,  £ = £^ , (£^+l) , . . . , (£^4 n) ) represented  by 
(2-1-6)  is  linearly  dependent. 
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It  is  concluded  that  the  complex  exponential  n-space  spanned  by  (2-1-1) 
is  invariant  under  the  differential  operator. 

The  waveform  (2-1-2)  is  now  operated  upon  by  the  integral  operator. 
First,  an  £-fold  reverse  time  integration  over  [t,°°)  is  performed  where 
t ^ 0.  Specifically,  it  is  readily  shown  that 
t t)  t. 


g(£) (t) 


^ )2  k m St  £ 

- .1. , - 


CO  CO  oo 


....  e n dt 

i=l  j-1  lj  q=l  q 


(2-1-7) 


k iik  A.  . 

y y1  ---A  -■ 

i=i  j=i 


£-1  j-1  SiT 
(t-x)  x e dx 


Ix(t) 


where  the  subscript  on  g denotes  an  £-fold  integration.  Concentrating 
on  I^(t)  and  making  use  of  the  binomial  expansion, 


1 p=0  p 


(2-1-8) 


I2(t) 


Integrating  12(f)  by  parts  successively  j+p  times  yields. 


£-1  £-1  » , j+p-1 

I (t)  = l (-1)P(  )t2~P_1  [ I (~l)qq!( 

p=0  P q=0 


j+p-1  j+p-q-1  St 


p=0 

£• 


.q+i 


e * ] 
(2-1-9) 


cl  J+E-1  p+q  4-1  j+p-1  tJ+*-q 

LJ  > Jq+r 


q— 2 Sxt 
e 


p=0  q-0 

Direct  integration  of  the  left-hand  side  of  (2-1-7)  reveals  that  it 
is  not  possible  for  (2-1-9)  to  contain  functions  in  t involving  powers 
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greater  than  (j-1).  With  reference  to  the  expression  in  (2-1-9),  this 
implies  for  integer  values  of  q less  than  (£-1)  that 


T (-l)p%!(t'1)(J+q'1)  - 0 

p=0  p q 


(2-1-10) 


The  net  result  is  that  (2-2-9)  can  be  rewritten  as 

1,(0  = l I (-l)P+qq !(V)(3+J  l)  e 1 

1 p=0  q=£-l  P * S?+1 


(2-1-11) 


Performing  a change  in  the  dummy  variable  q,  such  that  q = v+(£-l), 

£-1  j+p-£  , . j-v— 1 S.t 

Ix(t)  = I f (-I)?  (^-D!(  . )C'  ) Hi  6 \ (2-1-12) 

1 P=o;  =o  p ^_1  s ^ 


Finally,  substituting  (2-1-12)  into  (2-1-7)  yields 


J-v-1  S.t 


>(£) 


<t)  - l I (‘r1)  f e'1'. 

i=l  j-1  U l) • p=0  p v=0  v+*'  1 S ^ 

(2-1-13) 


From  (2-1-13),  it  is  concluded  that: 

(a)  Equation  (2-1-2)  is  infinitely  integrable  over  (°°,0]. 

(b)  g^(t)  is  a linear  combination  using  the  original  basis  in  (2-1-1) 
for  any  £ (i.e.,  space  invariance). 

(c)  For  any  integer  value  of  {g(£)(t)},£  = £^(£^+1) (^n-l) , 

spans  exactly  the  same  space  as  does  the  set  (2-1-1). 

(d)  Any  u+1  set  (i.e.,£  »£^,(£^+l,  ...,(£^+ii))  represented  by  (2-1-13) 
is  linearly  dependent. 

It  is  concluded  that  the  complex  exponential  n-space  spanned  by  (2-1-1) 

is  invariant  under  the  integral  operator  over  [0,°°)(i.e.,  the  results 
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for  operations  over  [0,°°)  and  (“,0]  are  related  through  the  sign  operator 

a 

(-1)  for  an  £-fold  integration) . 

The  limits  of  the  integral  operation  are  now  chosen  to  obtain  a 
finite  arbitrary  interval  of  integration.  Let  this  interval  be  [t  ,t]. 


g^(t)  in  (2-1-7)  becomes 


8(£)(t;to)  =J 


ti  \2  k m 


l t *„ 

J 1 i-1  i=l  1 

t t t J 


j-1  Sitl  J 
e II  dt 


q=l 


o o o 


(2-1-14) 


k m A.,  , St 

= l l 7ILTfr  f (t-T)£_1Td'  e 1 dr. 

i=l  j=l  ' ’’  J 


I3(t) 


Use  of  the  binomial  expansion  on  I^(t)  results  in 

1,(0  = (-l)P(£'1)t£“P_1[I9(t)  - I (t  )]. 

3 p=0  P 2 2o 


(2-1-15) 


It  follows  that 
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j+p-1 


1,(0  = l (-l)P(£“1)t£_P  1 [ l (-l)qq!(J+P_1) 
p=0  P q=0  q 


tj+p-q_leSit_  (j+p-q-leSito 


,q+i 


j+«.-q-2  sit 


p=0  q=0 

j+p-q-1  i o 
t e o i 

q-71 " 


,q+l 


(2-1-16) 


} 


Notice  that  the  finite  lltnlts  of  integration  have  introduced  into 
(2-1-16)  a polynomial  in  t of  order  £-1.  Again,  as  for  (2-1-9),  the 
coefficients  for  q<  £-1  are  identically  zero.  Hence,  following  the 
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steps  taken  in  the  development  of  ,2-1-13),  (2-1-14)  becomes 
k m.  A_  £-1  „ , j+|-£ 

v=0 


O-v-1  S . t tJ+p-^-VeSiCo 

t 1 r O ,.^-p-l-, 

•'  ^ 6 " [ IvSE ]t  } 


(2-1-17) 


,v+£ 


Examination  of  (2-1-17)  reveals: 

(a)  g,..(t;t  ) is  a linear  combination  of  the  functions  in  the  it+£  basis 

(£)  o 

S . t 

{t^  ^e  * ,l,...,t^  for  any  £ where  i = 1,2,... ,k;  j = 1,2,. ..,m^; 


l m = n. 
i=l 


(2-1-18) 


(b)  The  new  space  consists  of  the  subspace  spanned  by  (2-1-1)  plus 

£-1 

the  subspace  spanned  by  (l,...,t  }. 

We  next  consider  the  operation  of  passing  the  waveform  g(t),  as 
given  by  (2-1-2),  through  £ cascaded  identical  low-pass  filters.  Let 
the  impulse  response  of  each  filter  be  given  by 

At 

h(t)  = ■(  ’ (2-1-19) 


0 


t > 0 
t < 0 


where  A is  a negative  real  number.  The  impulse  response  of  the  £ cas- 
caded low-pass  filters  is 


ho  (t) 


/ £-1  At 
It  e 

(£-1)! 

0 


t >_  0 
t < 0 


(2-1-20) 


Passing  g(t)  through  the  cascade,  assuming  the  filtering  begins  at 


F 
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t = tQ,  results  la 


k An  f l-i  l-l  (S-TX)t 

«u>  <“'<>•»>  ' ,1  irHr ‘ I <”>  v 6 dT-  (2-l-21> 

-1  t 

o 


By  direct  comparison  of  (2-1-21)  with  (2-1-14),  it  follows  from 


(2-1-17)  that 


g(£)(t;to’X)  " ^ (£-1)!  pJ0  P J0  ( 1 (v  1 1)-(\H-S,-1) 


(S  -X)t 

_ j+p-£-v  i o 

i-V-1  S t tJ  p e w 0 „ , 

ft  i r o 1 _At._£  p 1 \ 

{ e - [— -SfZ ] e t } 

(S1-A)  (Si  - X) 


Notice  that  g^(t;tQ,x)  is  a linear  combination  of  the  functions  in 


//\  - _ » N 

,C  -i.  4.‘  ».  / 


the  iH-i.'  basis 


St  l-l  . 
{tJ_1e  1 , t 1 eAt}, 


i = 1,2, ... ,k 


j - 1,2,... ,m.  i=l 

= 1,2, ... ,£ 


[ mi=n 


(2-1-23) 


For  X = 0,  this  basis  reduces  to  that  of  (2-1-18). 


The  next  operation  to  be  con'  dered  is  that  of  passing  g(t)  as 
given  by  (2-1-2),  through  £ cascaded  identical  filters  having  single 
poles.  Let  the  linear  transfer  function  of  each  filter  be  given  by 


“r1 

n (S  - z ) 

P(S)  ,i=i 

Q(S)  = n 

n1  (s  - x.) 

j=i  J 


(2-1-24) 


It  is  assumed  that  the  real  part  of  each  pole  is  negative  and  that 
complex  poles  appear  in  conjugate  pairs.  The  transfer  function  of 
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the  £ cascaded  filters  is 


Hj,(S)  = [H(S)  ] 


ni~^  £ 

£ m j~i  J 

ni  j>" 

n-  (s  - x r 

j-i  J 


(2-1-25) 


Note  that  each  of  the  r>1  poles  of  fT(S)  is  of  £th-order.  Taking  the 
inverse  Laplace  transform  of  H^(S),  the  impulse  response  of  the  £ cas- 
caded filters  is 


V‘>  ■ j)  KS-x/  H((S>eSt]|s,Xr  (2-1-26) 


r 1 r £-1  i' 

I TJTJTT  I (j,  -l^r 

r-1  1,1  £ -1  \ 1 r r 


(£-£,)  £-1  X t 

. 1 r 
t e 


where 


Rr(s)  - (S  - \rr  H£(S) 


(i-ij) 


, th 


(2-1-27) 


and  (Xf)  denotes  the  (£-£^)  derivative  of  Rr(S)  evaluated 

for  S « Xf.  Passing  g(t)  through  the  cascade,  assuming  the  filtering 


begins  at  t • t , results  in 
o 


g(£)(t 


k mi  "l  A..  £ . (£-£,) 

(V 


(2-1-28) 


f £,-1  . . X (t-x)  SiT. 
j (t-T)  1 TJ_1e  r e dT  . 


To  carry  out  the  integration  in  (2-1-28)  we  may  use  the  results 
in  (2-1-22)  which  arose  from  (2-1-21)  since  the  kernels  of  both  Integrals 
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have  the  same  form.  The  result  is 


k mi  nl  A. 


CA-A, ) 


-*■  ~ n i V ~ — ~ i i 

!tf)(tito’^ = ^ TFI)T  a Rr  (V 


j 

J. 

*,-1  i -1  J+P“£i  P+\H-£  -1 

1(1  ) I (-D  (v+£  -l)!(^  ) (2-1-29) 

p=0  p v=0  1 1 


j+p-£  -v  (S  -A  )t 

j-v-1  St  t e r 0 A t £ -p-1 

( —6 — e 1 - U r t 1 


v+£. 


v+£ 


(S,  - X ) 
i r 


(Si  - Ar) 


} 


“jp^vri 


The  basis  representing  (2-1-29)  is  given  by 
, , St  £ -1  At  £ -1  Ant 

r . J “ J-  1±  1 1 1 i 

it  e ,t  e ,...,t  e } 


(2-1-30) 


where 


i=l,2,...,k  > k 

f l " 

j “ 1,2 mi  ) i=l 


1.2,. ...A 


The  new  space  formed  through  the  £-fold  operation  is  n+n^  dimensional 

and  includes  the  subspnce  (2-1-1).  If  we  let  t -*  °°  under  the  condi- 

o 

tion  Real  [S^  - A^]  < 0 for  all  possible  values  of  i and  r,  then  all 

of  the  coefficients  C . . are  identically  zero  and  the  £-fold  reverse 

jp£1vn 

time  operation  is  space  invariant. 


2-11 


Another  operator  that  preserves  space  Invariance  is  the  pure  delay 
operator.  If  we  operate  repeatedly  on  (2-1-2)  £ times  with  a delay 
operator  of  A seconds,  there  results 


g(t ;£A) 


k m, 

I 2* 


i*l  j-1 


1-1  i 

6(t-£A-T)xJ  e ax 

0 


km  . . -S  (t-£A) 

= l r A ,(t-£A)j  e 1 
i«l  j-1  J 


(2-1-31) 


k m.  j-1  • S.£A  , -S.t 

- 2 I*  A„  l 1 , 

i-1  j-1  J q=0  q 


t>£A>0  . 


Equation  (2-1-31)  verifies  the  space  invariance  under  the  repeated  delay 
operation. 


2 . 2 Correlation  Consider at  ions  Involving  Complex  Exponent i als 

In  general,  identification  of  the  complex  exponential  components 
of  a signal  involves  solution  of  a simultaneous  set  of  equations.  As 
the  correlation  between  the  components  is  increased,  the  equations 
become  more  ill  conditioned.  Consequently,  the  correlation  coefficient 
between  the  two  complex  exponentials  is  a measure  of  the  difficulty 
by  which  the  two  components  may  be  resolved.  By  way  of  illustration, 
two  signals  with  unit  correlation  result  in  a singular  set  of  equations 
whereas  two  s10nals  with  zero  correlation  yield  an  uncoupled  set. 

Consider  two  simple  complex  exponentials  given  by 
Sit 

(e  } ; i - 1,2,  0 < t < » (2-2-1) 

where  - 0^  +ja)j,  Sj  * - jou^,  < 0 and  j * /-l.  The  correlation 
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coefficient  between  exp  [S^t]  and  exp  [S^t]  over  the  time  Interval 
[T^,  VA]  Is  defined  to  be 

Tl+A 

Sit  V 

e L e dt 


P12(T1,T1+A)  T +A 


1 “ * 
s t St 

[ I e 1 e 1 dt 


T +A  * 1/2 

s2t  S,11 

e 1 e L dt] 


(2-2-2) 


Performing  the  integrations  and  rearranging  terms  yields 


P12(T1-V4) 


+ ao  + i (w,-u0) 


1 2 


e«,1«2)Sel(»1-»2)A_  1 J(4|  „ )T 

■_ , — - . .. ■ ...  e 

2T  A 2o„A 

(1  - e 1 )(1  - e 1 ) (2-2-3) 


In  general,  the  correlation  coefficient  is  a complex  quantity.  The 
initial  time  instant  T simply  adds  a constant  to  the  phase.  Since 
we  are  primarily  interested  in  the  magnitude,  there  is  no  loss  in 
generality  by  assuming  = 0.  For  this  special  case,  the  coefficient 
is  denoted  by  p^CA). 

Of  particular  interest  is  the  sinusoidal  case  for  which  = 0, 
i = 1,2.  Taking  the  limit  as  CT^  approaches  zero  with  T^  = 0,  (2-2-3) 


becomes 


P,,<A> 


3 A 
e - 1 


12  al,a2-0  j(^1-W2)A 


to.  - <o? 

sin( 5 A)  / 


(2-2-4) 


(0  - co„ 

1 2 


(O,  - (On 


The  phase  and  magnitude  of  the  correlation  coefficient  are  plotted 
in  Figure  (2-2-1)  where  it  is  assumed  that  w.  > Observe  that  the 


PHASE 


BOVNO 


Fig.  2-2-1.  Phase  and  magnitude  of  p (A) [ 


Hence,  for  an  observation  interval  of  length  A,  the  two  components 


are  uncorrelated  for 


with  n = 1,  we  obtain  the  convent icnal  condition  for  resolution  of  two 


frequency  components  which  is  given  by 


f,,  the  magnitude  of  the  correlation  coeffi 


For  other  values  of  f 


cient  is  bounded  by 
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P12(A) 


(2-2-8) 


°1?2*°  - 


1 * l2U 


This  bound  is  plotted  in  Fig.  2-2-1.  l/hen  , the  correlation 

coefficient  is  unity.  This  is  to  be  expected  since  the  two  components 
are  then  identical. 

The  situation  is  more  complicated  when  < 0;  i = 1,2.  This 
is  evident  by  consideration  of  the  spectra  involved.  For  the  infinite 


interval,  = 0 


yields  a line  spectrum  whereas  0.  < 0 results 


in  a continuous  spectrum.  Define 


a = 


3 = 


W1  " “2 


Y = A 


(2-2-9) 


where,  for  convenience,  it  is  assumed  _<  < 0.  Substitution  of 

(2-2-9)  into  (2-2-3),  with  T^  = 0,  the  co  relation  coefficient  is  ex- 


pressed as 


Pl2(A)  = 


/4a 


e(l+a)Yej6y  _ x 


1 + a + jB 


\J  [1  - e2Y 


(2-2-10) 


Ml  - e2aY] 


For  the  infinite  observation  interval,  A = 00  and  y = -°°.  This  yields 

fig ~ . (2-2-11) 


P12(oo)  1 + a + jB 


Observe  t.Knt 


pi2(A)  = -P12(~)  [ 


e(l+a)yejBy  _ x 


(2-2-12) 


\/[l  - e2Y][l  - e2aY] 


The  second  factor  in  (2-2-12)  can  be  interpreted  as  a "windowing"  factor 
due  to  the  finite  observation  interval. 


4 
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The  magnitude  and  phase  of  (2-2-10)  are  given  by 


~ o (l+a)Y  a 7 2 (l+a)Y 

1 - 2e  cos  gy  + e 


|Pl2<4)1'  VW  + 8 2 V 1 - (e2^+e2c"b  + e2  (1+^)  Y 


(2-2-13) 


/ Pl2(A)  - tan’1  [.(l+S)-  BY  _ 

L 8 sin  gy  + (l+a)(cos  8y  - e ^ ^ ) 


Because  this  is  a rather  complicated  expression,  attention  is  first 
devoted  to  the  case  of  the  infinite  observation  interval.  The  effect 
of  windowing  (i.e.,  finite  observation  interval)  is  discussed  later. 

As  is  evident,  either  from  (2-2-11)  or  (2-2-13),  the  magnitude  and  phase 
of  the  correlation  coefficient  for  A = 00  are 


|p12(“)|  = 


4a 


2 2 
(1+a)  + 8 


(2-2-14) 


/ / \ “1  r 0 

/P12(00)  = tan  I~~ 


For  a fixed  value  of  8,  the  maximum  value  of  | p (°°)  | occurs  for 


“m  = / 1 + e2  . 


The  peak  value  of  Ip^t00)!  is  then  given  by 


(2-2-15) 


I P12  (00>  1 1 c^ccn 


1 + 


A + 62 


(2-2-16) 


It  is  interesting  to  note  that  (2-2-10)  is  unchanged  by  interchanging 
the  subscripts  1 and  2 in  the  definitions  of  (2-2-9).  For  this  reason. 


For  |6|  » 1, 


lp12<")lla-l  “Tf-"  f to12<“5 1 la_a  1 


(2-2-18) 


n 


On  the  other  hand,  for  ( 8 1 <•■'  ], 


lp12(0o)lla»i  = 1 = 


(2-2-19) 


Also,  for  a >>  1 and  a >>  |b[,  observe  from  (2-2-14)  that 


lp12Mi=-  _ 


(2-2-20) 


Hence,  for  a fixed  value  of  | 6 | , the  correlation  coefficient  approaches 
zero  in  the  limit  as  a approaches  infinity. 

Equation  (2-2-14)  is  now  investigated  as  a function  of  8.  For  a 


fixed  value  of  a,  the  maximum  value  of  occurs  for 


8m“°- 


(2-2-21) 


It  is  interesting  to  note  that  the  correlation  coefficient  peaks  when 
(*>1  = but  does  not  necessarily  peak  when  = 0^  (i.e.,  a = 1).  For 


8„  = 0,  the  peak  value  of  the  correlation  coefficient  is 

n 
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lp12(oc,)  I *e=BM=0  l + a 


(2-2-22) 


For  a >>  1,  the  peak  becomes 


I p (oo)  I I S — 

i2  ^ ' 1 1 (3=8  =0 

M /a  • 


(2-2-23) 


Assymptotically , for  |B|  » a,  (2-2-14)  reduces  to 


|P12(“) 


(2-2-24) 


Therefore,  for  a fixed  value  of  a,  the  correlation  coefficient  approaches 
zero  in  the  limit  as  |S|  approaches  infinity.  By  comparison  of  (2-2-20) 
with  (2-2-24),  it  is  seen  that  the  correlation  coefficient  asymptotically 
approaches  zero  at  a faster  rate  with  respect  to  |B|  as  opposed  to  a. 

A plot  of  | 2 (°°)  I versus  |B|,with  a as  a parameter  is  shown  in 

Fig.  2-2-2.  Recall  that  the  larger  is  the  correlation  between  signal 
components,  the  more  ill-conditioned  are  the  equations  which  arise  in 
the  identification  problem.  Figure  2-2-2  points  out  that  the  problem 
of  resolution  is  eased  under  situations  of  either  small  a's  and  large 
B's  or  vice  versa.  Recognizing  that  a is  constrained  to  be  greater 
than  unity,  small  a implies  = a^.  Then  a large  value  of  |B|  is 
desirable  so  that  the  difference  in  and  u>^  will  aid  in  discriminating 
between  the  two  components.  On  the  other  hand,  if  |B|  is  small,  a 
large  value  of  a is  desirable.  It  is  interesting  to  note  that,  for  large 
a,  the  correlation  coefficient  is  relatively  insensitive  to  |8(  (e.g., 
see  curve  with  a = 100).  This  is  reasonable  since  a 1-^rge  value  of  a 
implies  that  one  component  decays  much  more  rapidly  than  the  other. 


(°°)  | as  a function  of  1 8 1 


Hence,  the  correlation  coefficient  is  more  influenced  by  the  relative 


decays  as  opposed  to  the  relative  oscillations.  In  the  identification 


problem  a and  6 are  specified  and  the  observation  interval  is  finite 


Since  an  infinite  observation  interval  was  assumed  in  obtaining  the 


curves  in  Fig.  2-2-2,  they  serve  as  a lower  bound  on  the  correlation 


coefficient  for  the  finite  interval 


A second  way  of  viewing  Ip^^C00)!  1®  presented  in  Fig.  2-2-3 
where  the  magnitude  is  plotted  as  a function  of  a with  | 6 | as  a parameter 


The  conclusions  arrived  at  from  the  previous  figure  are  still  valid 


(l-OI 
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The  contour  plot  is  shown  in  Fig.  2-2-4  where  the  parameter  is  |pi2(°°)I« 
This  plot  enables  the  user  to  determine  the  various  combinations  of  a 
and  IbI  which  result  in  a given  value  of  the  correlation  coefficient.  The 


plot  also  allows  one  to  determine  the  sensitivity  of  the  signal  parameters 
to  small  changes  in  the  correlation  coefficients.  For  example,  assume 


Fig.  2-2-4.  Equal  correlation  magnitude  contou’ s in  the 

a,  | 8 | plane.  a > 1 is  assumed  in  the  text. 
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0^  = -3  and  | @ | = 4.  For  | ^ (°°)  I = 0.4,  the  allowable  value  of  a 
is  approximately  23,  as  read  from  Fig.  2-2-4.  The  corresponding  value 
of  O ^ is  -69.  When  Ip^^00)!  = 0*5,  the  allowable  values  for  a are  1.3 
and  13.  This  yields  values  for  of  -3.9  and  -39,  respectively.  In 
this  case,  the  value  of  is  seen  to  be  highly  sensitive  to  changes 
in  the  correlation  coefficient. 

Our  discussion,  thus  far,  has  dealt  with  the  infinite 
interval.  We  now  consider  the  effect  of  a finite  observation.  The 
windowing  effect  is  accounted  for  in  (2-2-13)  by  the  second  factor 
which  is 


w (a , $ , Y ) 


tl  - 2ea+a)\os  BY  + e2U+a^ 
1 - (e2Y  + e2aY)  + e2(f+a)Y“ 


(2-2-26) 


w(a,B,y)  can  be  shown  to  approach  unity  for  all  possible  choices  of 
a and  8,  as  |y|  -*■  00 . Obviously,  the  effect  of  windowing  is  negligible 
when  w(a,B,Y)  - 1.  Since  a _>  1 and  y < 0,  v(a,8,Y)  is  bounded  by 

1 - e2Y  < w(a,B,y)  < (1  - 2e2Y)"1.  (2-2-27) 

Observe  that  the  bounds  are  independent  of  a and  3.  The  inequality  in 
(2-2-27)  can  be  used  to  obtain  an  estimate  for  the  minimum  length  of 
the  observation  interval  in  order  that  the  correlation  between  the 
two  components  be  approximately  the  same  as  for  the  infinite  interval. 

(In  general,  windowing  tends  to  increase  the  correlation).  From 
(2-2-27)  the  following  table  is  obtained: 
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Table  2-2-1. 

Lower  and  Upper 

Bounds  on  w(a,S,Y). 

Y 

l-e2Y 

(1  - 2e 

-1 

.864 

1.372 

-2 

.982 

1.038 

-3 

.997 

1.005 

t is  concluded  from  Table  2-2-1  that  the  record  length  can  be  assumed 
r.o  be  .finite,  as  far  as  w(C»,o,Y)  is  concerned,  provided  Y < -2.  In 
other  words,  when  the  length  of  the  observation  interval  is  such  that 

A > - — where  o£  < < 0;  1 = 2,3 (2-2-28) 

then  )p^(A)|  - Since  a finite  interval  tends  to  increase 

the  correlation  coefficient,  it  is  desirable  that  the  inequality  in 
(2-2-28)  be  satisfied. 

The  behavior  of|p,2(A)jas  a function  of  | Y I is  illustrated  it 
Figures  2-2-5  and  2-2-6.  In  Figure  2-2-5  a is  constrained  to  bi 
as  the  parameter  |ft|  is  varied  from  1 to  100.  Note  that  the  curve- 
have  settled  down  to  their  asymptotic  behavior  for  |y|  > 2.  It  is 
interesting  to  compare  the  damped  case  (i.e.,  i 0)  to  the  sinusoidal 
case  (i.e.,  » 0).  For  0^  - -1,  |S|  - | - u>2 1 . The  dashed  line 

in  Fig.  2-2-5  corresponds  to  the  sinusoidal  case  where  | » 10. 

This  serve.,  as  a reference  for  the  damped  case  where  ( g ( » 10  an.,  > ^ ■ 

-1.  In  Figure  2-2-6  |g|  is  constrained  to  equal  the  value  3 as  the 
parameter  a is  varied  from  1 to  100.  Again  asymptotic  values  have 
been  "reached"  for  |'  | > 2.  The  dashed  curve  in  Fig.  2-2-6  represents 


0.00  4.00  8.00  12.00  16.00  20.00  24.00 

tXICH) 

Fig.  2-2-5.  |pj,(A)  | vs . | Y | with  a = 1 and  |8|  as  a parameter. 
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v<^-  io  ^ 
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0-00  4.00  8.00  12.00  16.00  20.00  24  00 

(X  1 0-») 

Fig.  2-2-6.  |pi2(A)|  vs.|yI  wlth|8|  = 3 and  a as  a parameter. 
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the  sinusoidal  case  with  |co^  - | = 3.  Since  1 8 1 = 3,  each  curve 

in  the  figure  may  be  compared  with  the  sinusoidal  case  provided  is 
assumed  equal  to  -1.  For  a > 1,  it  is  interesting  to  note  that  there 
exist  values  of  |y|  where  the  correlation  coefficient  is  smaller  than 
that  of  the  sine0' idal  case.  The  only  exception  is  for  a = 1.  Notice 
that  the  largest  asymptote  occurs  for  a = J\  + 0 2 = /lO,  as  predicted 
by  (2-2-15). 

As  a final  topic,  the  orthogonality  of  sets  of  complex  exponen- 
tials is  discussed.  In  general,  two  simple  complex  exponentials, 
as  given  by  (2-2-1) , are  not  orthogonal  over  any  given  interval  with- 
in [0,°°).  Since  the  orthogonality  property  can  be  useful  in  certain 
situations,  it  is  desirable  to  be  able  to  orthogonalize  a set  of 
complex  exponentials.  The  Gram-Schmidt  orthonormalization  procedure 
is,  of  course,  applicable.  However,  a simpler  procedure  for  exponen- 
tial functions  was  developed  by  Kautz  [10]  and  applied  in  [11].  The 
orthonormalization  is  carried  out  over  the  semi-infinite  interval  and 
is  based  upon  the  Parseval  relation  for  an  inner  product  between  two 
time  functions.  Specifically,  the  Parseval  relation  is 


r 

fg*dt 

—00 


F (s)G*(-s) 
c 


_ds_ 

2*j  ' 


(2-2-29) 


It  follows  that  orthogonality  in  the  time  domain  is  equivalent  to 
orthogonality  in  the  frequency  domain.  If  f and  g are  sums  of  complex 
exponentials,  rhey  will  be  orthogonal  provided  F(s)G*(-s)  is  a rational 
function  which  is  analytic  (i.e.  , has  no  poles)  either  in  the  left 
half  or  right  half  of  the  s-plane.  Kautz  considered  the  set  of  expo- 


nent ials 


2- 


V 

(e  } , i = 1,2, .. . ,n t •>  0,  o < 0 (2-2-30) 

and  constructed  an  orthonormal  set  where  the  ntK  orthonormal  basis 
function  is  given  by 

n~1  * 

* 


r * , 11  (S  + V n-1  S + S* 

¥ (S)  * , -S  -S  ^ f*  n ---1- 

n '.nnn  S-S., S-S, 

n (s  - s > n 1=1  1 

i=l  1 


(2-2-31) 


To  verify  that  (2-2-31)  is  the  nC^  representative  of  an  orthonormal 
basis,  we  show  that 


vs>  y-s>  fj  ■ 


f i * = q 


i q 


(2-2-32) 


For  £ < q,  the  integral  in  (2-2-32)  becomes 


£-1 


1 * I 

V "VM  “VSq 

27Tj 


1ao  11  <S+Sk  > 11  11  (S-S.  )(-l) 

r k ,-i  i >.-r  2 k -£+i 


q-i 

r 

£_ 

£ q * 

-j°°  n (s-s  ) n (s+s  )(-i)q 

-lj-1  1 i2=l  2 


>q-1 

ds 


(2-2-33) 


Note  that  the  integrand  is  analytic  in  the  left-half  plane  and  on 
the  imaginary  axis.  Therefore,  closing  the  contour  to  the  left,  the 
resulting  integral  is  zero.  For  £ > q,  the  integral  is  analytic  in 
the  right-half  plane  and  on  the  imaginary  axis.  By  closing  the  contour 
to  the  right,  the  resulting  integral  is  again  zero.  For  £ = q, 

(2-2-33)  reduces  to 
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\j- S„-S 


rsi 


K- 


27Tj 


* 
S 

3-JL 


-1 


ds 


(-ix-W 


= i 


(2-2-34) 


S.  +S, 
k k 


-joo  (S-S4)(S  + S^) 

where  the  contour  can  be  closed  either  to  the  right  or  left  without 
affecting  the  result.  Equation  (2-2-31)  can  be  interpreted  in  terms 
of  passing  the  n^  exponential  time  function  through  an  all-pass 
filter  structured  from  the  previous  n-1  exponents.  The  all-pass  filter 
interpretation  points  out  that  it  is  the  phase  which  is  responsible 
for  the  orthogonality. 

Thus  far,  orthogonality  has  been  considered  over  a semi-infinite 
interval.  Provided  a finite  interval  is  suitably  long,  orthogonality 
can  still  be  approximated  by  the  Kautz  procedure.  From  (2-2-31),  note 
that 

, n S .t 

ip  (t)  = £ (S)}  = l A . e 1 . (2-2-35) 

n n . , l 


In  the  time  domain  orthonormality  requires 

m t 

1 . £ = q 

0 . £ 4 q . 


dt  = 


(2-2-36) 


Substitution  of  (2-2-35)  into  (2-2-36)  yields 


(;  (OiWOdt  = l l A A* 
q i-1  k-1 


r (s^t 

e dt 


£ t,  A A, 
V V 1 k 


L»  Lt  "h 

i-l  k-1  SJ  + S, 
i k 


1 

0 


£ - q 
£ 4 q 


(2-2-37) 


For  a finite  interval  of  length  A,  consider 


2-27 


ip0(t)>p  (t)dt 
l q 


f I *A 


i=l  k=l 


(S  « )A 
e - 1 

si  * < 


Clearly,  (2-2-38)  reduces  to  (2-2-37)  provided 


‘W4 


« 1. 


Assume  |a^|  < |o^|  for  i = 2,3,...,n.  Since 


(S  +S  )A 
i k 


< e 


2a^A 


(2-2-39)  can  be  replaced  by  the  similar  inequality 


2°1A  2v 

e = e ' « 1. 


(2-2-38) 


(2-2-39) 


(2-2-40) 


(2-2-41) 


Provided  the  interval  length  is  chosen  such  that  (2-2-41)  is  satisfied, 
the  orthonormal  basis  generated  by  (2-2-31)  is  very  close  to  being 
orthonormal  over  the  finite  interval.  Interestingly  enough,  |y|  ^ 2, 
which  was  the  condition  for  Ip^CA)!  = Ip^C00)!.  also  satisfies  the 
orthogona1 icy  condition  of  (2-2-41).  In  this  dissertation, 

M I 2 (2-2-42) 

is  the  condition  to  be  satisfied  if  a finite  interval  is  to  be  con- 
sidered as  though  it  were  an  infinite  interval. 
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2 . 3 A General  Parametric  Identification  Approach 

Consider  a linear  system  with  input  x(t)  and  output  y(t)  as 
shown  in  Fig.  2-3-1.  Assuming  the  system  is  described  by  a proper 
rational  transfer  function,  the  impulse  response  h(t)  may  be  expressed 
as 


h(t)  = 


f k m St 

i 1 

< i-i  i-i  J 

o 


t > ° 

t < o. 


(2-3-1) 


xlt) 


Fig.  2-3-1.  Linear  system. 

In  Section  1.3  it  was  show'  that  an  irreducible  representation  for 
this  system  is  the  input-output  relation  given  by  the  differential 
equation 


I D y(1)(t)  = l N x(l)(t).  (2-3-2) 

j-0  J i=0  1 

For  a complete  specification  of  the  equation  the  initial  conditions 
at  some  arbitrary  instant  of  time  tQ 

y '(t0);  k = 0,1,2 n-1 

x(L)(tQ);  Z-  0,1,2,.  ...m, 
must  also  be  specified. 


Notice  that  (2-3-2)  can  be  rewritten  as 


[ D y(j)(t)  - 1 N.x(l) (t)  = 0 . 

j=0  J 1=0 


(2-3-3) 


Equation  (2-3-3)  represents  a linear  combination  of  the  input  and  output 
signals  and  their  derivatives  up  to  the  orders  rr  and  n,  respectively. 


This  implies  that  the  set  of  functions 


(y(j),  x(i>; 


j = 0,1,2, ... ,n 
i = 0,1,2,... ,m 


(2-3-4) 


is  linearly  dependent.  Assuming  none  of  the  coefficients  in  (2-3-3) 
is  zero,  all  of  the  (n+m+2)  functions  in  (2-3-4)  are  needed  to  satisfy 
the  differential  equation.  Linear  independence  results  when  one  of 
the  functions  is  removed  from  the  set  (2-3-4)  since  it  is  no  longer 
possible  to  express  any  one  function  as  a linear  combination  of  the 
others.  This  suggests  that  (2-3-3)  is  defined  on  a space  having  dimen- 
sion (n+m+1). 

The  general  solution  of  the  differential  equation  given  by  (2-3-2) 

is  y^  ‘ (t)  = y(t)  and  consists  of  the  homogeneous  solution  y (t)  plus 

h 

the  particular  solution  y^(t).  The  impulse  response  h(t)  can  be 

determined  from  the  homogeneous  solution  alone.  Hence,  it  is  desirable 

that  any  identification  procedure  be  able  to  separate  y (t)  and  y (t) 

h p 

from  y(t)  conceptually,  if  not  in  practice. 

Since  a proper  rational  transfer  function  is  assumed,  m in 
(2-3-2)  is  strictly  smaller  than  n.  In  the  identification  problem  both 
m and  n are  unknown.  It  is  possible  to  simplify  determination  of  the 


values  for  m and  n by  first  assuming  that  m equals  n.  The  value  of  n 
is  then  found.  The  value  of  m follows  from  the  fact  that  the  coefficients 
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N.;  i = are  then  found  to  be  ~ero.  To  accomplish  this 


procedure  we  consider  the  set 


(v(J)  x(1)}  • J=l,2,....n 

’ i = 0,1,2 m n. 


(2-3-5) 


Note  that  yv  has  been  removed  from  (2-3-4),  as  mentioned  previously, 
this  results  in  linear  independence  of  the  set 


j = 1,2 n 

i = 0,1,2,. . . ,m. 


(2-3-6) 


Also,  observe  that  the  set  (2-3-5)  contains  the  function  ^,...,x^n^ 
which  were  not  originally  contained  in  (2-3-4).  The  set  (2-3-5)  will 
also  be  linearly  independent  provided  none  of  the  functions  >+1),..., 
x(n)  can  be  expressed  as  a linear  combination  of  the  functions  in 
(2-3-6).  The  functions  in  (2-3-5)  can  then  serve  as  a basis  to  span 
a (2n+l)  dimensional  space.  For  those  cases  where  (2-3-5)  constitutes 
a linearly  independent  set,  the  value  of  n can  be  determined  as  ex- 
plained next. 

For  purpose  of  illustration,  we  consider  an  idealizeo  identifi- 
cation problem.  Assume  for  all  t > 0,  that  analytical  representations 
are  available  for  the  system  input  and  output.  Consider  the  (2n+2) 


functions 


x(n)  (0) 

*•••**  9 J 9 


which  are  generated  as  shown  in  Fig.  2-3-2.  Assume  (2-3-5)  consti- 
tutes a linearly  independent  set.  The  object  is  to  determine  n and 
then  to  identify  the  system. 
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Fig.  2-3-2.  Reconstruction  oc  the  differential 
equation  describing  the  system  h. 


Because  known  analytical  expressions  are  assumed  for  x(t)  and 
y(t),  this  example  represents  a highly  simplified  situation.  Neverthe- 
less, it  serves  to  lliistrate  the  point.  The  procedure  begins  by 
assuming  the  system  differential  equation  is  given  by 


l D y(J)(t)  - l N.x(1)(t)  - 0 (2-3-7) 

j-0  J i-0 


Clearly,  (2-3-7)  remains  unchanged  when  multiplied  by  a constant. 
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For  convenience,  assume  the  equation  has  been  scaled  such  that  Dq  = 1. 
Equation  (2-3-7)  may  then  be  written  as 

r 

-y(0) (t)  = -y(t)  = l D.y(j)(t)  - l N.x(i)(t).  (2-3-8) 

j=l  J i-0  1 

The  functions  of  y^^(t)  and  x^(t)  are  known  for  all  values  of  j 
and  i.  For  a specified  value  of  r,  (2r+l)  of  these  functions  serve 
as  a basis  for  approximating  -y(t).  The  coefficients  can  be  deter- 
mined by  means  of  any  of  the  approximation  techniques  discussed  in 
Appendix  A. 

One  of  the  problems  is  to  determine  a suitable  value  for  r.  For 
r < n,  the  basis  is  inadequate  for  exactly  representing  -y(t)  and  a 
nonzero  approximation  error  is  involved.  When  r = n,  -y(t)  is  an 
exact  linear  combination  of  the  basis  and  the  approximation  error  i9 
zero.  To  demonstrate  how  the  order  of  the  system  is  determined,  we 
use  the  simplest  approximation  technique  in  Appendix  A which  is  that 
of  interpolation. 

Consider  the  (2ri-l)  equations  obtained  by  evaluating  (2-3-8)  at 
(2r+l)  different  time  instants.  This  results  in 

r r 

-y(0)(t  ) = 7 D y(j)(t  ) - 7 N x(1)(t. );  (2-3-9) 

j-1  j i“0  k 


k = 1 2r+l  . 


Equation  (2-3-9)  can  be  expressed  in  matrix  form  as 
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...y<r><tj) 

x(1) 

(1>(t2> 

<t2> 

"*(o)(t2> 
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y U2r+1) 


(t.) 
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(t2> 


• y (r) (t  x(0)(t  x(1)(t 

y 'c2r+l)  x '•c2r+l)  U2r+l)' 


.y<r)(tl) 

■x(r) (t2) 


•x(r)(t 


y(0)(tl> 

iy(0)(c2) 


|--- 

li"No 


!-N, 


2r+l)  -N 


; (t), 

P (t2r4-l)J 

« 

(2-3-10) 


The  square  matrix  is  nonsingular  provided  r < n.  Consequently,  n can 
be  determined  by  successively  increasing  the  value  of  r by  unity  until 
the  matrix  becomes  singular.  This  will  first  occur  when  r = n+1.  (The 
matrix  is  then  of  size  2n+3).  Having  determined  the  value  of  n,  (2-3-10) 
is  solved  for  the  coefficients  and  where  r is  chosen  equal  to  n. 
Under  the  assumption  that  (2-3-5)  is  a linearly  independent  set  the 
coefficients  for  i = (m+1) , ....  n will  be  found  to  be  identically 
zero  since  the  original  differential  equation  (2-3-2)  is  a linear  combi- 
nation of  the  basis  in  (2-3-6)  which  spans  a subspace  of  the  space 
spanned  by  (2-3-5). 

The  assumption  of  linear  independence  for  x^,  where  (w+1)  < i £ n, 


is  now  relaxed.  The  simultaneous  set  of  equations  in  (2-3-10)  are 
generated  as  before.  However,  in  addition  to  checking  the  singularity 
of  the  square  (2r+l)  x (2r+l)  matrix,  the  singularity  of  the  submatrix 
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(0)(tl) 

...  x(r)(t1) 

(0)(t2) 

xU)(t2) 

...  x(t)(t2) 

(0).  . (1),  . (r). 

x (tr+1)  x (tr+1)  ...  x (tr+1) 


is  also  checked.  As  long  as  the  submatrix  remains  nonsingular, 
while  the  value  of  r is  successively  increased  by  unity,  the 
procedure  remains  the  same  as  before.  The  procedure  is  changed 
should  the  submatrix  become  singular.  Suppose  this  occurs  for 
the  value  of  r equal  to  r^.  Then  both  the  (2r^+l)th  column 
and  row  are  removed  from  the  square  matrix  in  (2-3-10).  The 
reduced  matrix  is  checked  for  singularity.  If  this  matrix 
is  singular,  the  order  of  the  reduced  matrix  is  (2n+2)  and  the 
value  of  n has  been  found.  If  the  reduced  matrix  is  nonsingula  , 
the  procedure  proceeds  by  adding  to  the  reduced  matrix  an 
additional  column  to  the  [y^]  submatrix  while  adding  an 
additional  row  to  the  whole  matrix  as  shown  below. 


<ri> 

...  y (tj) 

(rj+l) 

y 

rr 

*-* 

!lx(0)(t  ) 

II  'V 

y(1)(t2) 

(ri) 

...  y 1 (t2) 

(rx+l) 

y 

(t2) 

!!."»(.,> 
il  1 

(rj+1) 


y (t5.  (t2r1+l)  y ^2^+1*  i 


2^+ir 


(r.-l) 

• x (tl) 

(r.-l) 

• x 1 (t,) 


(r.-l) 

•*  (t2r1+l) 


rl+l 


-N„ 


-N. 


-N 


v1) 


(0),  , 
y (tj) 


y(0)U2) 


y(0)(* 
y U2r;+1 
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This  process  is  continued  p times  until  the  resulting  (2r^+p)  x (2r^+p) 

matrix  is  singular.  The  order  of  this  matrix  is  r^  + n + 1 from  which 

n can  be  determined.  The  coefficients  a-e  found  by  solving  the  set 

of  equations  which  result  when  the  (r|+p)C*1  column  and  (2r^+p)t*1  row 

are  removed  from  (2-3-11).  This  procedure  also  yields  the  value  of 

m since  the  coefficients  N , , . ..,  N , will  be  identically  zero. 

nr*-!  r^-i 

The  previous  discussion  utilized  a differential  equation  descrip- 
tion of  the  system.  This  was  not  necessary.  Many  other  representa- 
tions exist  for  a dynamical  system.  Provided  each  representation  is 
equivalent,  any  given  input  will  result  in  the  same  output  irrespec- 
tive of  the  representation  used.  Equivalency  between  the  convolution 
integral  and  differential  equation  representations  was  demonstrated 
in  Chapter  1. 

A variety  of  equivalent  representations  can  be  generated  from 
(2-3-3)  by  operating  on  the  differential  equation  with  various  linear 
operators.  Assume  a particular  linear  operator  is  denoted  by  g^. 

Applying  this  operator  to  (2-3-3),  there  results 

gj  l D y(J)(t)  - l N x(1)(t)}  = 0 . (2-3-12) 

j-0  J i=0  1 

Making  use  of  linearity,  as  defined  in  (1-1-1),  (2-3-12)  can  be 
written  as 

I D [g  y(J)(t)]  - l N [g  x(i)(t)]  - 0 . (2-3-13) 

j-0  J i-0  1 k 

This  procedure  can  be  extended  to  a cascade  of  £ such  operators. 
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Let  g^(t)  denote  the  impulse  response  obtained  when  the  unit  impulse 
is  applied  to  the  operator  g . The  Laplace  transform  of  g (t)  is 

K K 

denoted  Dy  0 (S) . Since  convolution  in  the  time  domain  corresponds 
to  multiplication  in  the  frequency  domain,  the  process  of  operating 
on  (2-3-3)  by  the  sequence  of  operators  g2>  g^  can  le 

described  by 


t-t 

n r o 


„ t-t 
m ( o 


where 


I Dj  J f j^(T)y''  ’ (t-x)dx  - l ] fc(T)xU,(t-T)dT  - 0 ; t > t_  (2-3-1A) 


(2-3-15) 


t»(t)  - i."1  I n G (S)J. 


For  convenience,  define 


v<(t;*) 


Xjfr  n 


f^(T)y(j)(t-T)dT 


f j (T)x(1) (t-T)dT. 


(2-3-16) 


Equation  (2-3-14)  can  be  rewritten  then  as 

n m 

l D y (t;fc)  - In.  x . (t ; «.)  - 0 
j-0  J J i"0 


(2-3-17) 


where  it  is  understood  that  the  system  is  initially  relaxed  at  time 
t • The  representation  in  (2-3-17)  is  in  terms  of  the  new  set  of 
functions 
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x^t ;£)) 


j = 0,1,... ,n 

i = 0,1 m. 


(2-3-18) 


Provided  these  functions  span  an  n + m + 1 dimensional  space,  the 
representation  in  (2-3-17)  is  equivalent  to  that  in  (2-3-7).  Note 
that  the  coefficients  and  are  identical  in  the  two  representa- 


tions. 

Equation  (2-3-17)  has  significant  implications  as  far  as  the 
identification  problem  is  concerned.  In  general,  only  finite  records 
of  the  input  and  output  are  available.  The  differentiation  scheme 
shown  in  Fig.  2-3-2  for  generating  a set  of  n + m + 2 functions  to  be 
used  in  the  identification  is,  in  general,  impractical.  Differentia- 
tion tends  to  enhance  the  noise  and  measurement  errors  inherent  in 
experimental  data.  Alternative  schemes  are  suggested  by  (2-3-17) 
which  incorporate  the  n + m + 2 functions  of  (2-3-18),  In  other 
words,  operators  other  than  differentiators  can  be  utilized  to  generate 
a suitable  set  of  n + m + 2 functions  from  the  input-output  data. 

A practical  difficulty  yet  remains  with  use  of  (2-3-17).  The 
functions  in  (2-3-16)  are  generated  from  derivatives  of  the  input  and 
output.  In  practice,  only  the  input,  x^°\  and  output,  y^°\  are 
given.  The  question  then  arises,  "Is  it  possible  to  generate  an 
equation  similar  to  (2-3-17),  where  the  coefficients  and  are 
preserved,  but  the  functions  involved  are  generated  only  from  the 
given  input  and  output?" 

To  answer  this  quesfion,  consider  the  frequency-domain  input- 


output  relation  for  a nonrelaxed,  linear,  time-invariant,  causal  system 
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which  is  given  by  (1-3-6).  Recall  that  the  coefficients  C arise  from 

q 

the  nonzero  initial  conditions  at  time  t = 0.  Rewriting  (1-3-6),  there 
results 


n m n-1 

l D SJY(S)  - l NS  X<S)  - l CSq  = 0 (2-3-19) 

j=0  J i-0  q-0  q 


Let  the  transfer  function  of  a linear  cascade  of  Z filters  be  denoted 
by 

Z 


F (S)  = l n Gk(S)]. 
v k=l 


(2-3-20) 


Multiplication  of  (2-3-19)  by  F^(S)  yields 


n . m , n-1 

l D SJF  (S)Y(S)  - l NS  F (S)X(S)  - [ C SqF.(S)  = 0.  (2-3-21) 
j=0  J i=0  1 £ q-0  q £ 


The  inverse  Laplace  transform  of  F^  (S)  , as  pointed  out  in  (2-3-15), 
is  f- (t).  Hence, 

p 

£ _1{SPF^(S) } = ^ {f4(t)J  - f£p)(t).  (2-3-22) 

dt 

Converting  (2-3-21)  to  the  time  domain,  we  obtain 


l D1  I y(t)f(j)(t-T)dx  - l N I 

j=0  J J £ i-0  1 J 


x(T)f^1) (t-T)dT 

0 


T C fjq)(t)  - 0 
q-0  q 


(2-3-23) 


where  it  is  understood  that  differentiation  is  with  respect  to  the 
variable  t.  For  convenience,  define  the  functions 
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Yj (t - y(T)qJ'(t-T)dr 


(2-3-24) 


xA(t;i)  - j x(T) f ^ (t-x)dt  . 

0 

Notice  the  reversed  roles  of  the  functions  in  (2-3-16)  and  (2-3-24). 
Equation  (2-3-23)  now  becomes 


1 D y (t;£)  - [ N x (t;£)  - £ C > (t) 


(2-3-25) 


For  the  relaxed  case  the  coefficients  C are  identically  zero  and 

q 

(2-3-25)  is  identical  in  form  to  (2-3-17).  For  the  system  representa- 
tion in  (2-3-25)  to  be  equi.alent  to  the  or.m.,.  'inferential  equa- 
tion* ^(t)  should  be  'n  differentiable  and  t tie  •-••• 


(fjp;(t)}  ; p = 0,1 n 


(2-3-26) 


should  span  an  (n+l)-dimensional  function  space.  A block  diagram  of 
(2-3-25)  is  given  in  Figure  2-3-3. 

In  general,  arbitrary  inputs  and  nonzero  initial  conditions 
must  be  considered.  For  simplicity,  a relaxed  system  whose  input  is 
a unit  impulse  is  first  discussed.  For  a system  with  impulse  response 
given  by  (2-3-1),  the  output  is 


(km.  , S.t 

l ‘ 

l-U-l  ‘J 

0 ; t < 0 


(2-3-27) 


Noti  :e  that  y(t)  is  a linear  combination  of  the  n basis  functions  in 
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Input  Output  Initial  Conditions 


I 


0 

Fig.  2-3-3.  Block  diagram  of  (2-3-25). 


i.2-1-1) . Because  the  output  is  the  impulse  response,  we  are  dealing 
with  the  homogeneous  case.  Hence,  it  is  necessary  to  treat  only  the 


output  signal.  Any  space-invariant  linear  operators  which  generate 
a suitable  function  f^  (t)  may  be  used.  The  coefficients  can  be 
determined  by  using  one  of  the  techniques  in  Appendix  A to  construct 
from  the  (r+l)-basis  set 


{yi (t;Z)} 


j - 0,  . . . , r 
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a set  of  (r+i)  homogeneous  linear  equations  In  the  unknowns  . The 
next  step  Is  to  determine  the  value  of  n,  as  discussed  earlier  in 
this  section.  The  linear  equations  can  then  be  reduced  to  n non- 
homogeneous  equations  by  selecting  Dq  equal  to  unity.  The  coefficients 
Dj  are  used  to  write  the  characteristic  equation  of  the  system  from 
which  the  poles  are  extracted.  The  resulting  n complex  exponential 
time  functions  are  then  used  to  fit  the  output  record  y(t)  to  deter- 
mine the  residues  in  (2-3-27).  Knowledge  of  the  residues  and  poles 

is  sufficient  to  determine  the  coefficients  N^. 

Previous  techniques  described  in  the  literature  can  be  inter- 
preted in  this  context.  For  example,  the  differential  operator  is 
used  in  [23],  [24],  This  procedure  is  usually  applied  to  idealized 
synthesis  problems  and  is  impractical  for  the  identification  case  where 
measured  records  are  involved.  References  [25],  [26]  make  use  of  the 
time-reverse  integral  operator  given  in  (2-1-7).  For  finite  records 
truncation  errors  are  incurred.  This  error  can  be  reduced  by  using 
a suitably  long  interval,  as  shown  in  Section  2.2. 

A nonrelaxed  system  with  a general  input  x(t)  is  now  considered. 

The  time  response  is  the  inverse  Laplace  transform  of  (1-3-6).  This 
results  in 

km.  . . S.t  km.  , . S.T 

y(t)  = l l B tj_1e  1 + x (t— T ) [ l l A TJ  e ]dx  (2-: 

i=l  j-1  i-1  j-1  1J 


where  the  coefficients  8. , are  functions  of  the  coefficients  C , which 

lj  q 

are,  in  turn,  functions  of  the  initial  conditions  and  the  coefficients 


-28) 


Dj  and  N^.  From  the  point  of  view  of  (2-3-i9),  it  is  of  interest  to 
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identify  the  2n  + m + 1 parameters 


{D  :N  ;C  } 

j i q 


j = 1... 
i = 0,.. 
q = 0... 


• » n 

• , m 

. , n-1 


(2-3-29) 


The  response  in  (2-3-28)  can  be  interpreted  as  the  sum  of  two  responses. 
Let 


' k m-  1-1  V 

e 


hi(t)  = 


i ii 


i-i  j-i 
0 


t > 0 


t < 0. 


(2-3-30) 


Hence,  one  of  the  responses  is  the  response  of  a relaxed  system  with 
impulse  response  h^(t)  to  a unit  impulse.  The  second  response  is  the 
response  of  a relax'd  system  with  impulse  response  h(t),  as  given  in 
(2-3-1),  to  the  input  x(t).  This  interpretation  is  illustrated  in 
Fig,  2-3-4.  Notice  that  both  h^(t)  and  h(t)  utilize  the  same  basis 
functions  given  in  (2-1-1). 


J(t)( 


k 


| V1 


X (t)o- 


k It) 


c+ 

J 


'Jit) 


Fig.  2-3-4.  Interpretation  of  (2-3-28). 


2-4  3 


Identification  of  the  system  requires  generation  of  the  set  of 
functions 

:(q) 


(f^qy(t);  y . (t ;£) ; x.(t;£)}  ; j = 0,  ....  r 

1 1 i = 0 r. 


(2-3- 


The  coefficients  in  (2-3-29)  can  be  determined  by  using  one  of  the 
techniques  in  Appendix  A to  construct  from  the  (3r+3)-basis  set  of 
(2-3-31)  a set  of  (3r  + 3)  homogeneous  linear  equations  in  the  unknown 
coefficients.  The  next  step  is  to  determine  the  value  of  n,  as  dis- 
cussed earlier  in  this  section.  It  is  unnecessary  to  check  the  func- 
t • ' !it)  for  linear  independence  because  they  are  chosen  to  be 

linear lv  independent,  as  explained  in  conjunction  with  (2-3-26).  In 
general,  the  coefficients  are  not  of  interest  because  they  depend 
on  initial  conditions  which  are  input  dependent.  If  desirable,  the 
(3r  + 3)  equations  can  be  reduced  to  a set  of  (2r  + 2)  equations  by 
means  of  a systematic  elimination.  The  remaining  equations  can  then 
be  solved  for  the  coefficients  D.  and  . 

Thus  far,  the  identification  approach  has  utilized  operators 
which  preserve  the  numerator  and  denominator  polynomial  coefficients 
of  the  linear  syste  ' rational  transfer  function  (i.e.,  the  coefficients 
and  Nj  of  the  differential  equation).  This  is  not  necessary.  Opera 
tors  which  do  not  preserve  these  coefficients  can  be  used  proviied  a 
transformation  can  be  determined  which  enables  the  original  differen- 
tial equation  coefficients  to  be  obtained  from  the  new  coefficients. 
This  results  in  a much  larger  class  of  operators  than  those  previously 


31) 


discussed. 
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In  analogy  with  Fig.  2-3-3,  consider  the  block  diagram  of 
Fig.  2-3-5.  The  impulse  response  of  the  filters  in  the  filter  banks  of 
Fig.  2-3-5  arc  no  longer  required  to  be  successive  derivatives  of  each 
other  as  was  the  case  in  Fig.  2-3-3.  In  fact,  it  is  not  even  required 
that  f^(t);  k = 0,  ,.,,n  be  differentiable.  The  coefficients 


<V  V V ! 


j 0 , • • • , n 
i = 0,  . . . , m 
q = 0 , . . . , n-1 


(2-3-32) 


are  chosen  to  satisfy  the  constraint 


n m A n-j.  „ 

l 6 F (S)Y(S)  - V N F (S)Y(S)  - T C F (S)  = 0 

j=0  j J i=0  1 1 q=0  q q 


(2-3-33) 


where  F^(S)  is  the  Laplace  transform  of  f^(t).  The  system  identifica- 
tion can  be  based  upon  (2-3-33)  provided  a transformation  can  be 
developed  between  the  "hatted"  coefficients  of  (2-3-32)  and  the  "unhatted" 
coefficients  of  (2-3-29).  Of  course,  the  block  diagram  of  Fig.  2-3-5 
reduces  to  that  in  Fig.  2-3-3  when  f^(t)  = f ^ it). 

For  this  purpose,  consider  the  mapping  of  variables 


Fk(S)  = Vk  where  S = F^k(vk) 


(2-3-34) 


and  a unique  inverse  is  assumed  to  exist.  (Multiple-valued  inverses 
are  permitted  provided  a meaningful  unique  inverse  car.  be  defined.) 


Use  of  (2-3-34)  in  (2-3-33)  results  in 


? - " n v^if'1^1)]  - "y1 

j=0  J J i-0  q=0 


C = 0.  (2-3-34) 
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Finally,  define 


Y(v)  = Y[G_1(v) ] and  X(v)  = X[G-1(v)]. 


Equation  (2-3-36)  can  then  be  written  as 


n A ^ n-1 

l D.V1?^)  - l N V X(v)  - l 
j =0  J i=0  1 q=0 


C V — 0 . 

q 


(2-3 


(2-3 


X(v)  and  Y(V)  can  be  interpreted  as  the  input  and  output,  respectively, 
of  a system  whose  transfer  function  in  the  v-plane  is 


H(v)  = 


l N V* 
i=0 


(2-3 


l V 


j 


j=0 


Having  determined  H(v),  the  poles  and  zeros  of  H(S)  can  be  obtained 

A 

from  the  poles  and  zeros  of  H(v)  through  the  transformation 


S = G_1(v). 


(2-3 


Knowledge  of  the  poles  and  zeros  of  H(S)  is,  of  course,  sufficient 
to  construct  the  numerator  and  denominator  polynomials  and,  therefore, 
the  coefficients  and  . 

As  a final  point,  the  mapping  of  (2-3-35)  suggests  that  the 
block  diagram  realization  of  (2-3-33)  is  more  efficiently  realized 
through  the  cascade  connections  shown  in  Fig.  2-3-6. 


-37). 


-38) 


-39) 


-40) 
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2 . 4 Comparison  of  the  General  Parametric  Identification 
Approach  with  Different  Operators 

As  discussed  in  Section  2.3,  a variety  of  operators,  cur  be  used 
in  the  general  parametric  approach  to  system  i d 'nt  i f Icat ion.  Since 
the  system  impulse  response  is  assumed  t.  b*  ,i  1 near  combination  of 
complex  exponentials,  as  in  (2-11  , • • ' the  property 

that  an  exponential  subspace  is  pr  < ref  rated  operations  of 

linear  operators.  Under  idealize:  i is.  . results  obtained 

should  be  independent  of  the  operat  r • . practice  the  choice  of 

an  operator  can  affect  the  solutions  because  different  operators  result 
in  different  sensitivity  to  noise,  different  degrees  of  ill-conditioning 
of  the  matrices  involved,  etc.  Therefore,  a suitable  operator  should  be 
chosen  for  a given  situation.  Also,  a proper  approximation  scheme  from 
Appendix  A should  be  chosen.  Although  a trade-off  between  "quality"  of 
the  identification  and  ease  of  implementation  of  the  approximation 
scheme  usually  exists,  in  general,  an  increase  in  computational  complexity 
does  not  guarantee  an  increase  in  quality. 

We  start  by  first  introducing  a noiseless  ideal  example  where 
the  operations  are  performed  analytically.  The  system  is  assumed  to  be 
nonrelaxed  (i.e.,  an  unknown  excitation  is  assumed  to  exist  before  the 
identification  procedure  begins).  In  the  initial  p-oblem  the  input  is 
assumed  to  be  exponential.  Later,  an  input  of  a non-exponential  type 
is  employed.  Also,  noise  is  introduced  at  which  time  the  identifica- 
tion is  carried  out  entirely  numerically. 
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Example  2-4-1 

The  system  to  be  identified  is  shown  in  Fig.  2-4-1.  Let  its 
impulse  response  be  given  by 


2e_t  cos  t;  t > 0 l e(“1+i)t  + e("1_i)t  ; t > 0 


h(t)  = 


(2-4-1) 


; t < 0 


; t < 0 


The  probing  signal  for  the  identification  is  chosen  to  b< 
-2t 


x(t)  = 


; t > o 

0 ; t < 0 


(2-4-2) 


If  the  system  had  been  relaxed  at  t “ 0,  1 1 . « 


e C (cos  t + sin  t)  - e : t 


yx(t)  = 


; t • 0 


f-|  [(l-j)e(-1+j)t  + (1+j  )e  (" 
0 


However,  the  system  is  assumed  to  be  n<-  -relaxed  at  t \s-  unit  its 

response  due  to  the  non-zero  initial  conditions  resul  s in 


yx(t) 


e 1 (cos  t - sin  t)  ; t _>  0 
? ; t < 0 

f-|[(l+j)«("1+j),:  + (l“j)e(_1~^)t] ; t > 0 

: t < 0 


(2-4-4) 


Hence,  the  overall  output  as  a result  of  the  response  to  x(t)  and  the 
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initial  conditions  is  given  by 


Fig.  2-4-1.  A non-relaxed  system  to  be  identified. 

Several  different  operators  are  used  to  solve  this  example.  In 
each  case  the  input  and  output  are  given  as  in  (2-4-2)  and  (2-4-5) , 
respectively.  With  the  information  provided,  the  problem  is  to  iden- 
tify the  system. 

Fjr  purpose  of  illustration,  we  begin  by  taking  advantage  of 
the  fact  that  x(t)  is  exponential.  Hence,  it  can  be  included  within 
a modified  system  which  is  excited  by  an  impulse  as  presented  in 
Fig.  2-4-2.  The  object  now  is  to  determine  the  modified  impulse 


response  hm(t)  from  which  h(t)  can  be  developed.  The  solution  is 


Fig.  2-4-2.  The  modified  svs'em  h (t). 

tn 


Case  1 - The  original  Prony  method 

First  we  evaluate  the  signal  at  several  points  as  given  in 
Table  2-4-1.  This  implies  that  a sampled  version  of  the  given  data 
is  used.  The  sampling  increment  has  been  chosen  properly  to  omit 
the  occurrence  of  multiple  Riemann  surfaces.  The  sampled  version  of 
the  genera}  linear  combination  of  exponentials,  as  given  in  (2-1-2), 
is 
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q 

yq  (T  = .5) 

0 

1 

l 

.6966820193 

2 

.2621969375 

3 

-.01821986209 

4 

-.1309543389 

5 

-.1382616922 

6 

-.1010564004 

7 

-.05746896525 

8 

-.02427926367 

9 

-.004806869416 

10 

.003777201613 

Table  2-4-1 


With  reference  to  (2-3-33)  and  (2-3-35),  a cascade  of  shift  operators 
in  the  sampled  data  is  equivalent  to  the  mapping 


G(S)  = eST  and  Fk(S)  = [eST]k  = zk  (2-4-8) 

(the  variable  z is  used  in  place  of  V,  to  conform  with  the  usual 
notation  for  this  particular  mapping.)  The  equivalent  system  of 
Fig.  2-4-2  is  then  characterized  by  the  homogeneous  equation 


l D eSTJY(S)  = 0.  (2-4-9) 

J-0  J 


However , 

i-  _1{eST^Y(S) } = y (t  + jT) 


This  results  in  the  time-domain  equation 


n 

l 6 y(t  + JT)  = o. 

j-0  J 


(2-4-10) 
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Let  y(jT)  - y^  . For  t-0,  the  above  equation  becomes 


l K y,  - °- 


(2-4-11) 


Equation  (2-4-11)  may  be  rewritten  as 


I Yj  - - 

i=i  J J 


y where  D = 1. 


(2-4-12) 


It  is  interesting  to  note  that 
n-1  A 

yn“-  * Yj 

j=0  J J 


where  D = 1 
n 


(2-4-13) 


serves  as  an  extrapolation  schema. 

The  first  step  is  to  determine  the  value  of  n.  To  start  we  set 
r = 1 and  check  the  singularity  of 


1 

o 

yl 

1 

.6966820193 

yl 

y 2 

6966820193 

."621969375 

(2-4-14) 


The  computed  value  ot  • e determinant  is 


det.  Mj  - .2231688985. 


(2-4-15) 


Since  is  nonsingular,  we  set  r * 2.  The  corresponding  matrix  is 


yo 

yl\ 

y2  " 

1 

.6966820193 

.2621969375 

*1 

y 2 1 

y3 

- 

.6966820193 

.2621969375 

“.01821986209 

72 

* 1 

y3 

y4 

.2621969375 

".01821986209 

“.1309543389 

Checking  the  singularity  of  , it  is  found  that 


det.  M2  = .004211290373. 


(2-4-16) 


(2-4-17) 
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is  also  non-singular.  Hence  is  checked  for  singularity.  This 
time 


det.  = det. 


y0  yi'  ?2\  y3 


I I 


yl  y2 I y3'  y4 


y2  y3  y4  I y5 


y3  y4  y5  y6 


0 . 


(2-4-18) 


Having  determined  that  is  singular,  the  order  of  the  system  is 


established  at  n = 3.  We  now  set  equal  to  unity  and  solve  for  D^, 


D2>  using  the  following  set  of  equations: 


“ /\  “ 

D1 

V 

°2 

= - 

yl 

S . 

.y2. 

(2-4-19) 


The  solution  to  (2-4-19)  is 


A 

D1 

"5.612059902 

A 

°2 

= 

10.58438618 

A 



D. 

7.389056099 

3 

(2-4-20) 


Knowledge  of  the  coefficients  is  adequate  to  obtain  the  characteristic 
equation  of  the  difference  equation  in  (2-4-11).  In  particular,  the 
characteristic  equation  is 

1 /V  1 

0 . (2-4-21) 


A /N  /N  9 * 3 

D0  + D1Z  + D2z  + °3Z 
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The 


roots  of  the  characteristic  equation  in  the  z-phase  are 


Z1 

.3678794412  + j.O 

Z2 

= 

.5322807302  + j. 2907862882 

Z3 

.5322807302  - j. 2907862882 

(2-4-22) 


Using  the  inverse  of  the  mapping  in  (2-4-8) , we  have 


S = G_1(z)  = | <Ln  z- 


Hence,  the  poles  of  the  system  in  the  S-plane  are 


(2-4-23) 


V 

Z, 

1 

-2  + jO 

S„ 

= ^ £n 

z _ 

= 

-1  + jl 

2 

T 

2 

S3 

,Z3 

-1  - jl 

(2-4-24) 


The  natural  frequencies  in  (2-4-24)  art  identical  to  those  in 
y(t)  = h (t)  given  by  (2-4-5).  Observe  that  (2n+l)  = 7 samples  were 

IR 

needed  (i.e.,  y~  through  y ) in  order  to  establish  the  order  of  the 
0 6 

system  at  n = 3.  On  the  other  hand,  2n  = 6 samples  were  needed  (i.e., 

A A A 

y^  through  y^)  in  order  to  determine  the  coefficients  D^,  ^rom 

which  the  system  poles  were  computed. 

In  connection  with  Prony's  method,  two  questions  arise:  (1)  What 
should  be  the  sampling  rate?,  (2)  How  long  a data  record  is  needed? 

The  following  observations  are  pertinent: 

a)  Closely  spaced  samples  give  rise  to  highly  dependent  sets 


of  linear  equations.  A measure  of  the  dependence  is  the  magnitude  of 
the  determinant  of  the  set  of  equations.  For  this  example,  the 
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magnitudes  of  the  determinants  corresponding  to  r = 1 and  2 are  plotted 

in  Fig.  2-4-3  as  a function  of  the  sampling  interval  T.  Note  the  small 

magnitudes  which  occur  for  small  values  of  T . This  is  reasonable 

because  the  variation  in  the  signal  is  small  for  closely  spaced  samples 

2 

A column  dependence  occurs  even  if  one  is  willing  to  use  (n+1) 
different  samples  in  the  form 


y0 

yi 

y2 

...  V 

' n 

\ 

V1 

ynL+2 

. . . yn^+n 

5 nn  > • • • > 1*2  > 

s 

\+l 

yn2+2 

yn2+n 

yn  +1 

yn  +2 

yn  +n 

n 

n 

n 

n 

b)  Widely  spaced  samples  also  give  rise  to  highly  ill-conditioned 
sets  of  linear  equations.  Note  the  small  magnitudes  which  occur  for 
large  values  of  T in  Fig.  2-4-3.  This  arises  when  dealing  with  transient 
waveforms  due  to  the  large  dynamic  range  associated  wi th  the  initial 

and  final  samples  of  the  signal.  Hence,  the  final  rows  and  columns  are 
numerically  small  with  the  result  that  the  magnitude  of  the  determinant 
is  also  small. 

c)  An  optimum  sampling  interval  exists.  For  example,  the 
curves  in  Fig.  2-4-3  peak  for  T = 0.7.  By  maximiz*  < det.  M^|  with 


respect  to  T,  an  estimate  of  the  optimum  value  of  the  sampling  interval 
can  be  obtained.  Knowledge  of  this  value  can  also  be  used  to  obtain 
a rough  estimate  of  the  system  order  n.  From  Section  2.2  it  is  known 


fX  1 0-2) 


r 


Fig.  2-4-3.  Magnitude  of  the  determinants  corresponding  to 

r ■ 1 and  2 as  a function  of  the  sampling  interval 
T.  The  scale  for  |det  M | is  40  times  that  for 
|det.  | > 
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that  the  length  of  the  observation  interval, A,  should  be  chosen  such 
that 

Y = Act1  <_  -2.  (2-4-25) 

In  general,  for  a known  value  of  n,  2n  samples  are  needed  to  determine 
the  system  poles.  This  requires  an  observation  interval  of  length 
A = (2n  - 1)T.  Substituting  for  A in  (2-4-25),  there  results 

n-I["aT+1]-  (2-4-26) 

The  inequality  in  (2-4-26)  can  be  used  to  obtain  n,  a first  guess  for 
the  value  of  n.  This  is  done  in  the  following  way.  The  value  of  T 
is  chosen  to  be  that  which  maximizes  |det.  | . The  value  of  equals 
the  largest  real  part  of  the  system  poles  (recall  < 0)  and  is  not 
known.  However,  for  transient  waveforms,  the  over-all  rate  of  decay 
is  governed  by  the  slowest  significant  decaying  component.  Hence,  a 
guess  of  can  be  obtained  from  the  rate  of  decay  of  the  response 
envelope.  For  our  example,  with  T = 0.7  and  = -1,  n is  estimated 
to  equal  2.  This  is  close  to  the  true  value  of  3.  Having  made  an 
estimate  of  n,  it  is  >ossible  to  estimate  a suitable  value  for  A. 

Since  (2n+l)  samples  are  needed  to  determine  system  order,  an  interval 
of  length  2nT  should  be  used.  Therefore,  a suitable  value  of  A is 

A A 

A = 2nT. 


For  our  example,  A = 2.8. 


(2-4-27) 
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Case  2 - The  Overdetermlned  Prony  Method 

The  basic  difference  between  the  overdetermined  Prony  method 
[16],  [17],  [18]  and  the  original  Prony  method  [13]  is  in  the  type  of 
approximation  scheme  utilised.  While  the  former  utilizes  a pseudo- 
inverse technique,  which  is  equivalent  to  the  least-square  technique, 
the  latter  uses  an  interpolation  scheme.  (See  Appendix  A).  Let  n^ 
be  greater  than  n,  the  order  of  the  system.  With  reference  to  (2-4-11), 
an  overdetermined  set  of  equations  is  given  by 


’ y0  yl  y2  •••  yn 
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yn  yn  +1  yn  +2  yn  +n 

.0  . 

111  1 

(2-4-28) 


where  (n^  + n + 1)  different  samples  are  used.  To  solve  for  the  D^'s, 
we  premultiply  (2-4-28)  by  the  transpose  of  the  rectangular  matrix. 
This  yields 
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(2-4-29) 


Equation  (2-4-29)  can  be  rewritten  as 
n,  n 

2 1 ^l+ic^2,+i^?.  = ® » k = 0,1,2, 

i=0  <1=0  K * 1 K 


,n. 


(2-4-30) 


With  Dq  arbitrarily  set  to  unity,  (2-4-30)  can  be  rewritten  as 


n -1  n A n -1 

I l v1+ky£+1D£  - - xl  yi+ky±  k = 0,1,2 n-i.  (2-4-31) 

i=0  <1=1  1 * 1 1=0  1 K 


This  is  of  the  same  form  as  the  least-squares  equation  (A-l-7).  As 
with  the  original  Prony  method,  the  square  matrix  in  (2-4-29)  has  a very 
small  determinant  for  closely  spaced  samples  (i.e.,  a highly  linearly 
dependent  set  of  equations  results).  For  instance,  a record  of  y(t) 
of  length  5 with  T = .05  (i.e.  , 101  samples)  yields  the  following 
determinant  for  r = 2: 
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2 

98 

I y 

i-0 

98 
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Iwi 
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^ y i+ly i+2 
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= det . 


11.49921692 

10.97923768 

10.42420972 


10.97923768 

10.4992277 

9.984006238 


10.42420972 

9.984006238 

9.508731645 


= 2.029682836  • 10 


(2-4-32) 


This  Is  an  extremely  small  value.  In  general,  oversampling  results  in 
small  determinant  values. 

Oversampling  does  not  offer  any  advantage  even  with  more  widely 
spaced  samples.  In  the  original  Prony  method,  previously  discussed, 
it  was  pointed  out  that  a sampling  interval  of  T = 0.7  is  desirable. 

For  a modest  oversampling  ratio  of  R = n^/(r  + 1)  = 5,  n^  = 15  when 
r = 2.  The  corresponding  record  length  is  A = 2n^T  = 21.  As  pointed  out 
in  Section  2.2,  a record  length  of  2 is  adequate  for  = -1,  Hence, 
the  record  length  of  5 used  in  conjunction  with  (2-4-32)  already  repre- 
sented a very  long  record.  The  record  length  of  21  offers  negligibly 
more  information  than  was  contained  in  the  record  length  of  5 due  to 

21°i  Q 

the  exponential  nature  of  the  signal  (i.e.,  e <<  e <<  e ).  Hence, 
negligible  improvement  results  from  overdetermining  the  problem  in 
this  way.  Since  the  lower  rows  in  the  rectangular  matrix  of  (2-4-28) 


ire  extremely  small,  they  contribute  negligibly  to  the  inner  products 
<:  1-29) . Thus,  the  inner  products  obtained  for  the  oversampled 
it.  approximately  the  same  as  those  obtained  by  using  the  least- 
um  ■ technique  with  the  square  matrix  of  the  original  Prony  method. 
• ' r , the  usual  overdetermined  Prony  does  not  offer  any  advantage 
. original  Prony  method. 

Bv  introducing  a delay,  it  is  possible  to  achieve  an  advantage 
! v e r the  original  Prony  method.  This  newly  proposed  approach  is  dis- 
cussed next . 


Case  3 - The  Delay  Operator  Method 

Once  again,  Aet  n^  be  greater  than  the  system  order  n.  Consider 
the  overdetermined  set  of  equations  in  (2-4-33). 

fy„  y,  y,,  ...  y„  If6„l  To  1 


yo 

yJ 
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Xnj+1 
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yj+2 

y2J+2 

ynj+2 

\ 

yj-l 

>2j+n1 

ynj+n 

(2-4-33) 


Observe  that  the  original  data  record  must  be  of  length  A = (n,  + nj)T, 

1 d 

where  T^  denotes  the  sampling  interval.  Each  column  of  the  rectangular 
matrix  is  composed  of  (n^  + 1)  coi.aecutive  time  samples.  However,  the 
significant  difference  from  (2-4-28)  is  that  consecutive  columns  in 
(2-4-33)  are  advanced  by  j time  samples.  When  j = 1,  (2-4-33)  reduces 
to  (2-4-28).  Each  column  can  be  interpreted  as  a basis  function.  The 
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matrix  equations  can  then  be  viewed  as  a linear  combination  of  the 
basis  functions.  Whereas  the  basis  functions  in  (2-4-28)  do  not  change 
very  much  from  one  to  the  next,  they  are  drastically  different  in 
(2-4-33).  In  this  approach,  is  chosen  to  be  as  small  as  desired. 

The  value  of  j is  then  determined  from  the  relation 

jTd  = T (2-4-34) 


where  T is  the  preferred  sampling  interval  of  the  original  Prony  method. 
Wote  that  n^  = (A/T^)  - nj . 

A 

As  before,  the  D^'s  are  obtained  by  premultiplying  (2-4-33)  by 

A 

the  transpose  of  the  rectangular  matrix.  With  D = 1,  we  obtain  the  n 

ji  w^A  ■ - Wi  ■ k - 


We  now  illustrate  the  advantage  of  this  approach.  As  in  cases  1 and 

2,  let  the  record  length  be  A = 5.  Let  the  sampling  interval  T,  equal 

d 

.05  (this  was  also  the  sampling  interval  used  in  Case  2).  Hence,  101 
samples  are  available  as  ir.  "ase  2.  From  Case  1 it  is  known  that  the 
preferred  value  of  T is  0. 7.  Hence,  j = 14  and  n^  = 58.  For  r = 2,  the 
following  determinant  is  obtained  when  checking  for  the  system  order: 
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"1.116650923 

3.418188405  1.731065907 

.1945935105 

"1.116650923  .1945935105 

.5143890218 

(2-3-36) 


.06549016746 


From  Fig.  2-4-3,  the  value  of  the  determinant  for  the  original 


Prony  method  (r  = 2)  is  calculated  to  be  approximately  .0039.  From 

(2-4-32),  the  value  for  the  overdetermined  Prony  method  (r  = 2)  is 

appcoximately  .000002.  These  values  are  one  order  of  magnitude  and 

four  orders  of  magnitude  smaller,  respectively,  than  the  value  in 

(z-4-36).  A definite  improvement  is  noted  for  this  case  without  having 

expended  any  more  computational  effort  than  in  Case  2.  A plot  of  the 

continuous  basis  functions  (i.e.,  y(t),  y(t  + 0.7),  y(t  + 1.4),  and 

y(t  + 2.10)  is  given  in  Fig.  2-4-4  for  t > 0.  Note  that  these  basis 

functions  correspond  to  successive  columns  of  the  square  matrix  in 

(2-4-33)  where  y,  ...  = y(t  + (kj+i)T.)  is  evaluated  at  t = 0.  These 
kj+i  d 

basis  functions  involve  the  advance  operator  exp  [ST]  of  (2-4-8). 
Solution  of  (2-4-35)  for  the  coefficients  yields 

fD.l  ["7.135606018] 


16.54686248 


(2-4-37) 


16 . 4446467  7 


The  characteristic  equation  of  the  corresponding  difference  equation  is 


a a ~ 2 /v  3 

D0  + + d2z  + d3z  = o. 


Solving  for  the  roots  in  the  z-plane,  performing  the  transformation 
to  the  S-plane  results  in 


51  |zl]  “2  + J° 

;2  ■ (jTd)_1in  z2  - -1  + jl 

S z,  -i  " Jl  . 


(2-4-38) 


Once  again,  the  poles  of  H(S)  have  been  found. 
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Fig.  2-4-4.  Continuous  basis  functions  structured  through  the  delay  operator. 
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Prony's  method,  and  variations  thereof,  deal  with  samples  of  the 
observed  waveform.  One  could  just  as  well  work  with  the  continuous 
time  response,  as  appears  in  (2-4-10).  In  addition,  it  is  not  necessary 
to  utilize  advances  (l.e.,  negative  delays)  to  generate  a set  of  basis 
functions.  For  variety,  consider  the  delay  operator 

-ST 

G(S)  - e . (2-4-39) 


The  corresponding  time-domain  equation, 
n 


l D y(t  - jT)  = 0. 

j=0  3 


analogous  to 


(2-4-10), 


is 


(2-V40) 


The  set  of  basis  functions  generated  from  the  delay  iperator  is  also 
shown  in  Fig.  2-4-4  where  it  is  understood  that  y(t  - jT)  equals  zero 
for  t < jT.  With  reference  to  the  integrated-squared-error  approach 
in  Appendix  A (see  (A-l-11)),  the  n equations  for  the  coefficients  are 
given  by 


n 

I [ 

i=l 


y(t-kT)y(t-iT)dt]D 


y (t-kT)y (t)dt ; k 


iiT 


nT 


0*  1 ^ * • 
(2-4-41) 


>n-l 


where  the  k basis  function  is 
gk(t)  = y (t  - kT). 


(2-4-42) 


Since  the  largest  delay  is  nT  and  because  we  do  not  want  to  use  the 
zero  portion  of  a basis  function  in  the  inner  products,  the  integration 
is  over  the  inteval  (nT,A). 
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Case  4 - The  Reverse  Time  Integral  Operator 

The  parameters  of  h (t)  are  now  identified  using  the  reverse 
m 

time  integral  operator  introduced  in  (2-1-7)  to  generate  a set  of  basis 
functions.  This  operator  was  first  utilized  by  Ccrr  [25]  and  then  by 
Jain  [26]  who  developed  their  approaches  from  different  points  of 
view.  The  results  of  this  dissertation,  however,  make  clear  that  the 
two  approaches  are  essentially  the  same. 

The  simplest  development  of  the  identification  scheme  with  the 
reverse  time  integral  operator  proceeds  from  the  time-domain  version 
of  (2-3-21).  Let  the  operator  f {•}  consist  of  an  n-fold  reverse 
time  integration  from  00  to  t.  In  particular, 


1 tJ. 


OO  00  00 

1 


f4{,}  * / / *'•••/  Z {•>  q?idtq 


(2-A-A3) 


Carrying  out  an  4-fold  integration,  where  42:n,  the  time-domain  version 
of  (2-3-21)  is  then  given  by 


V D.  / y^(t)dr  - J N-  f 

1 1 ‘ (4.-1) ! y U;dT  1 "l  > (4-1) 

OO  OO 

j-0 


4-1 


x(i)(x)dr 


i-0 


(2-4-44) 


n-1 


- J c,  ^ -<I=T>T—  «<q)(^)<3T  - 0 


q-0 


where  6^^(t)  is  the  singularity  function  which  is  the  derivative  of 


the  unit  impulse,  6 ( t ) . Since  the  singularity  functions  occur  at  the 
origin,  they  are  not  within  the  interval  of  integration  and  each  term 
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in  the  summation  over  the  index  q is  identically  zero.  Let  the  i.-fold 
integrations  of  y^(t)  and  x^'(t)  from  « to  t be  deonted  by  y^  ^ (t) 
and  x^i-!l\t),respectively.  Equation  (2—4—44)  can  then  be  written  as 

n m 

l D y(J'A)(t)  - l x(i"°(t)  = 0.  (2-4-45) 

j=0  i=0 

Making  a change  of  variables  involving  the  indices,  (2-4-45)  becomes 

n m 

V D .y(n_4'j)(t)  - I N x(m_£"i)(t)  = 0.  (2-4-46) 

“ n-j  m-i 

j=0  i=0 

For  £=n,  (2-4-46)  yields 

n m 

l Dn_jy(_j)(t)  - l Nm_ix(m"n“i)(t)  = 0.  (2-4-47) 

J-0  i-0 

Hence,  under  the  reverse  time  integral  operation  the  coefficients  of 
the  differential  equation  description  of  the  system  are  preserved  in 
reversed  order.  Notice  that,  for  a case  where  only  exponential  signals 
are  involved,  one  could  arrive  at  the  same  conclusion  directly  from 
properties  (c)  and  (d)  of  (2-1-13). 

For  variety,  the  generalized  integrated-squared-error  approach  is 
utilized  as  the  approximation  scheme  for  generating  the  simultaneous 
equations  involving  the  coefficients  (see  Appendix  A,  (A-l-34)).  The 
values  of  L and  K are  chosen  to  be  2 and  1 respectively  (l.e.,  the 
coefficients  are  computed  such  as  to  minimize  the  sum  of  integrated- 
squarcd-errors  between  the  signal  and  its  approximation  and  the  first 
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time  reverse  integral  of  the  signal  and  its  approximation).  Assuming 
a system  of  order  r,  we  may  write 


0 r t 

1 l [/  2y^)(t)yj(e)(t)dt  ] Dr_.=0;  k-0,1,2,. 


. ,r.  (2-4-48) 


*=-l  j=0  L1 


In  matrix  notation,  (2-4-48)  yields  for  r=l 
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and 
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(«), 


(t)  - y(/'k)(t)  = 


k-l-l 


(k-i-i) : 


yQ0)(Odt 


The  function  y^^(t)  is,  of  course,  identical  to  y(t)  of  (2-4-5). 


For  r=l,  the  basis  set  needed  is  given  by 

(y(0)<O.  y^V). 


(2-4-50) 


Hence,  applying  a time  reverse  integral  operator  successively  to 
y(t)  in  v2-4-5)  yields  the  basis  functions 
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(2-4-51) 
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y(  2)(t)  = / y(-1)(T)dT  = l 


( -t  , -2t 

-e  sint  - l/4e  ; t^O 


; t<0 


Let  the  approximation  interval  [T^,  ] be  [0,5].  Substituting 


(2-4-51)  into  (2-4-49)  yields 

fr 


M, 


D. 


.1126700914  -.03122914547 

[-.03122914547  .1905598834 


.5499617419  -.1251736779 


-.1251736779  .1126700914 

(2-4-52) 


Computing  the  deteminant  of  results  in 


det.  M_ 


.1764679909. 


(2-4-53) 


Hence,  is  nonsingular.  For  r=2,  it  is  determined  from  (2-4-48)  that 


the  basis  needed  is  given  by 


<y(0\  y(-2).  ,M)). 


(2-4-54) 


The  new  basis  function  is 
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For  r=2,  the  determinant  of  M,  is 
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1 

t 

.1126700914 

-.03122914547 

-.03453295368 

= det.  < 

-.03122914547 

.1905598834 

-.1952578315 

-.03453295368 

-.1952578315 

.241469618 

’ .5499617419 

-.1251736779 

-.2374113928 

+ 

-.1251736779 

.1126700914 

-.03122914547 

-.2374113928 

-.03122914547 

.1905598834 

Hence,  M 

2 is  nonsingular. 

(-4) 

For  r=3,  y ’ 

has  to  be  added 

(2-4-56) 


=.0005575168764. 


basis  set.  Applying  the  time  reverse  operator  to  y 


(-3) 


yields 


.(-4), 


y ' ' (t)  = / yx  (T)dt= 


(-3) 


t cost  - ^ e "t:;  t>0 


t<0. 


(2-4-57) 


The  determinant  of  M is  computed  to  be  det.  M =0.  Hence,  the  system  order 

i 3 

is  established  at  n = 3.  The  basis  functions  utilized  for  the  GISE  approxim- 
ation are  shown  in  Fig.  2-4-5  over  the  time  interval  [0,5]. 

It  is  of  interest  to  compare  the  magnitudes  of  the  determinants  which 
arise  in  the  three  separate  cases  of  1)  an  ISE  approximation  to  y(t),  2)  an 
ISE  approximation  to  y^  ^(t),  and  3)  a GISE  approximation  to  both  y(t) 
and  y^  ^(t).  In  the  context  of  the  GISE  formulation,  these  cases  arise 
with  1)  L=l,  K=0,  2)  L=l,  K=l,  and  3)  L=2,  K=l,  respectively.  The  magnitudes 
of  the  determinants  for  r=l  and  2 are  tabulated  in  Table  2-4-2  for  these 
three  cases.  Note  that  case  3 results  in  larger  values  of  the  determinants 
for  the  nonsingular  matrices  and  is,  therefore,  superior  to  cases  1 and  2. 


(t-OIX) 
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MICROCOPY  RESOLUTION  TEST  CHART 


Table  2-4-2 


~ — ■ ■-  — ■ ■ i 

GISE 



r 

case  1 
L-l,  K-0 

case  2 

L-l,  K-l 

. . 

case  3 
L-2,  K-l 

1 

|det.  M.J 

.0463 

.0205 

.1765 

. 

2 

Idet.  M2I 

.0000791 

.00000494 

.000558 
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Both  Carr's  method  [25]  and  Jain's  method  [26]  employ  the 
reverse  time  integrator  to  generate  the  basis  functions  and  ISE  approx- 
imation to  generate  the  simultaneous  equations  for  the  coefficients. 

However,  different  solution  methods  are  used  on  the  equations. 

Consequently,  the  poles  converge  to  their  true  values  along  different 
trajectories  as  the  value  of  r is  increased  from  1 to  3.  By  way 
of  illustration,  the  equations  resulting  for  Case  3 (i.e.,  L-2, 

K=l)  are  solved  using  the  two  solution  methods  o Carr  and  Jain. 

Equation  (2-4-48)  may  be  rewritten  as 

l l [J5  ya)(t)yW(t)dt]  D = - f /5  ya)(t)ya)(t)ya)(t)dt  (2-4-58) 
l—l  j=l  0 k j r 3 i— 1 0 k k 0 

; k-0 , 2 , . . . , r-1 

where  has  arbitrarily  been  set  equal  to  unity.  With  Carr  i method,  (2-4-58) 
is  solved  directly  for  the  coefficients.  The  results  for  r'  2,  3 are  shown 
in  Table  2-4-3. 


Table  2-4-3 


B 

— 

D 

r 

' 1 - — - 

°r-l 

°r-2 

■2S 

m 

1 

4.236699946 

wm 

1 

1.496951865 

1.633885061 

3 

1 

4 

6 

4 

■■■ 

■■■■■■ 

■■■■■■■■■■ 
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Solution  of  the  corresponding  equations  for  r=l,2,3,  respectively  yields  the 
poles  presented  in  Table  2-4-4. 


Table  2-4—4 


r 

— ‘ " 1 1 — 1 H 

real 

imaginary 

1 

-4.236699946 

0 

2 

-.7484759325 

-.7484759325 

1.036179926 

-1.036179926 

3 

-1 

1 

-1 

-1 

-2 

0 

In  Jain's  method  the  r+1  coefficients  for  the  characteristic  equations  are 
given  by  the  square  roots  of  the  r+1  diagonal  co-factors  of  M^.  In  particular. 


(2-4-59) 


where  is  the  determinant  of  after  removing  the  f*ow  and  column. 
For  a specified  value  of  r,  the  characteristic  equation  in  terms  of  the 
diagonal  co-factors  is  given  by 


(2-4-60) 


The  resulting  coefficients  for  r=l,2,3  are  found  in  Table  2-4-5. 


Table  2-4-5 


r 

Dr  " 

Vi  ■ FZ 

°r-2  = J^22 

°r-3  “ \/A33 

1 

.5506632136 

.8140220103 

2 

. 2823259536 

.460784953 

.4200809338 

3 

.00590294882 

.02361179528 

.03541769292 

.02361179528 

Solution  of  the  corresponding  characteristic  equations  for  r-1,2,3  respectively 
yields  the  poles  presented  in  Table  2-4-6. 
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Case  5 - The  Simple  Integral  Operator 

A representation  of  h(t)  through  h (t)  is  possible,  in  the  identi- 

m 

fication  process,  only  for  exponential  inputs.  In  general,  it  is  not 
required,  and  in  this  case  the  identification  will  be  carried  out  directly 
on  the  system  as  represented  in  Fig.  2-4-1.  This  will  serve  as  an  intro- 
duction to  the  case  of  arbitrary  (i.e.,  nonexponential)  inputs. 

The  integral  operator  was  utilized  by  [31]  and  introduced  in  (2-1-14). 
In  particular 


{•}  = Jc  fh 

t t 

o 


J*2  {•}  ft 

to  q=l 


dt 

q 


= I C dx. 

to  (£-D! 


Substitution  of  (2-4-61)  in  (2-3-23)  for  £>n  yields 


(2-4-61) 


1-i-l 


* -i, -i 


n-l 

" I c 

q-0  < 


H-q-1 


x(t)  dx 

(2-4-62) 


(A-q-1) I 


dx 


From  (2-4-62)  it  is  understood  that  the  simple  integral  operator  preserves 
the  coefficients.  Let 


y'J  ~'(t) 


Jt  (£-j- 


■il 


itrli 


1)1 

£-1-1 


(£-1-1) 1 


y(t)  dT 


x(t)  dx 


(2-4-63) 


. . £-q-l  t-t  £-q-l  . (t  t0) 

t jt-xj 1_  dT  = J 0 U du  = 


£-q 


o a-q-Di 


(£**q— l)  l 


(*-q)  1 


Substituting  (2-4-63)  into  (2-4-62)  yields 


” D. 

j-0  J i=0  1 


:<1_l)(t)  - l C 


-1  _ (t-t0) 


£-q 


q=0  q (£-q)  1 


=0  . (2-4-64) 


With  a change  of  variables  involving  the  indices,  (2-4-64)  becomes 


l Dn-j  y 

j=0  n J 


m 


(n-£-j ) . y N x(o-£-i) 

4 In  0-1 


(t-t0) 


£+q+l-n 


(t) 


i-0 


n-1  

(t)  Cn-q-l  (£+q+l— n) ! 

q-0 

(2-4-65) 


For  £ = n and  t “0,  (2-4-65)  yields 

0 


n . . N m 

I D_  , y("j)(t)  - [ N 

j-0  n j i-0  m 1 


(m-n-i) 


n-1 

(t)  - I c 


v i 

b n"q"1 


tq+1 

(q+1)! 


0.  (2-4-66) 


Equation  (2-4-66)  suggests  that  the  coefficients  of  the  differential  equation 
are  preserved  in  reversed  order.  While  the  block  diagram  representing 
(2-4-62)  is  given  in  Fig.  2-3-3,  the  block  diagram  representing  (2-4-66)  is 


given  in  Fig.  2-3-6,  where  D - D , and  N - N , and  C , = C . 

n— r r tn— r r n— r— 1 r 
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From  (2-4-66) , the  basis  functions  involved  are 


{y'Wi<t),  ...,y''~n;(t),  x(m_n)(t),  X(m-n  1)(t),..,x("n)(t), 


(m-n-1) 


1!’  2!’  **’  n!  ; . 


(2-4-67) 


For  r = n = 2 and  m = 1,  (2-4-67)  becomes 

{y(0) (t) , y('1)(t),  y(-2)(t),  x(-1)(t),  x(_2)(t),  t,  —■  ) 

Then,  the  basis  functions  are  computed  to  be 


2-4-68) 


y(0)(t)  = y(t)  = 


— f _9  f- 

2e  Cos  t - e 


; t > 0 
; t < 0 


- / y<0)(r)dT  - ( 

f 1 -t,_ 

2 ~ e (C°s t - 

0 

l ? 

(1.1  -t 

('2>(0  ■ I - ! 

f 2 t + 4 " e Sint 

i ? 

; i i -2t 

- J dT- 

1 2 - 2e 

\ 

0 

(0 

1 -2t 

2 e 


x(_2)(t)  = / x(_1)(T)dT 
0 


lit  - ) + i.-** 


2 4 4 

0 


; t > 0 
; t < 0 
; t > 0 

; t < 0 
; t >.  0 
; t < 0 

; t > 0 
; t < 0 


along  with  the  polynomial  components  t and 


(2-4-69) 


The  ISE  approximation  scheme  yields  the  following  set  of  linear 
equations: 


<y(0\ 

0 

V 

<y(0). 

y(-1}> 

<y<0>. 

(-2)  1 (0) 
y > <y  , 

,„(0>  (~2)  1 (0)  , „ 
y f * 9 py  » c > 

< y(0\ 

t2. 
2 ‘ 

, y(0)> 

, y(-x)> 

<y^_1). 

y(-2)>|<y(-1), 

, x(-»> 

x<‘2V  t > 

< y^11. 

t2. 
’ 2 ' 

<y<-2>, 

, y(0)> 

<y(-2\ 

, y('1)> 

<y(-2>. 

y(-2)>|<y(-2>, 

, x(-«> 

<y('2).  x<^>  </-»■  t > 

< y(_2), 

t2 

’ 2 

77 

7«> 

77 

7-7 

77 

y<-»» 

,7-7 

Xx^.  x<-2>>|<x(-1>.  t > 

"77 

y(0)> 

<x<-2>. 

y(-l)> 

<x(-2). 

y(-2)>l<x(-2)> 

1 

*<-”> 

<x(-2>.x<-2)>(<x(-2).  t> 

<x(-2). 

<t.  y(0>  > <t.  /-« 


4 ■ y<°»  4 . yM)  > <V  • y(‘2)> 


4.*i'i)>  4"('2) 


>|  < t,  t> 

*\  4- ** 


.2  .2 

<v . v> 


(2-4-70) 


Equation  (2-4-70)  can  be  solved  directly  by  Jain's  method.  Or,  by  Carr's 
method  after  setting  D2  * 1,  moving  the  first  column  to  the  right  hand  side 
of  the  equation,  and  removing  the  last  row. 

In  general  this  is  an  insensitive  operator.  As  was  pointed  out  in 
section  2.1,  and  as  is  evident  from  (2-4-69),  an  additional  . '•lynomial 
space  is  structured  via  the  simple  integral  operator.  Since  polynomials 
are  ever-growing  functions,  even  over  short  records  they  become  the 
dominant  contributors  to  the  set  of  basis  functions.  This  contribution 
is  further  magnified  if  the  input  is  a decaying  exponential.  Thus,  the 
main  contribution  of  the  basis  functions  of  (2-4-69)  to  the  inner  products 
of  (2-4-70)  will  be  obtained  from  the  polynomial  part  of  the  basis  functions. 


2-80 
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Since  these  particular  basis  functions  include  only  a constant  and  a 
linear  term,  one  should  expect  that  the  set  (2-4-70)  will  be  highly  ill- 
conditioned.  Hence,  as  in  Case  2,  the  simple  integral  operator  is  another 
operator  that  generates  a highly  linearly  dependent  set  of  basis  functions. 

Next  we  consider  an  example  where  a nonexponential  input  is  applied 
to  the  system  and  additive  noise  is  present  at  the  output. 


Example  2-4-2 


The  system  to  be  identified  is  shown  in  Fig.  2-4-7.  Let  its  impulse 


response  be  given  by 
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From  (2-4-71),  the  system’s  transfer  function  is 


H(S)  - 


[s-(-2Vj)] 


1 . 
2j 


[S-('^-j)]  [S-(-~  + *5j)]  [S-(-i-.5J)] 


=2  • 


< W+2-1,’*  f.'> 

<S  +i  - j)(S+i  + j)  <S  + i - . 5J)  <S  + i + .51) 


2.  (S  + .4411)  (S  + .3478  - j 1.0086)  (S  + .3478  + jl.0086) 

(s  +h~ j)  (s  +h  + j)  (S  + ¥ " -5j)  (S  + i + *5J) 


(2-4-72) 


The  probing  signal  for  the  identification  is  chosen  to  be 


f sin  t 
x(t)  -i  C 

lo 


; t > 0 
; t < 0 


(2-4-73) 


The  system  is  assumed  to  be  relaxed  at  t“0.  The  output 


y(t)  - / x(T)h(t  - T)dT 

0 


(2-4-74) 


was  calculated  numerically  Over  the  time  interval  [0,6tt]  using  601  samples 
spaced  .Olir  seconds  apart  for  each  of  x(t)  and  h(t).  Figure  2-4-8  presents 
the  input,  the  output,  and  the  system's  impulse  response  over  the  time  in- 
terval [0,6it]. 
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The  noise  n(t)  is  pseudo-ranck ..  and  uniformly  distributed  between 
-A  and  A.  A noise  sample  function  of  101  samples  over  the  interval  [0,6tt] 
is  plotted  in  Fig.  2-4-9.  The  value  of  A is  chosen  such  that  the  signal  to 
noise  ratio  for  the  noise  corrupted  output  is  given  by 


101  „ 
l y\ 

10  1o810  -T5T-7 


= lod£ 


(2-4-75) 


where  101  samples  of  each  of  y(t)  and  n(t)  were  considered  over  the  interval 
[0,6tt].  The  samples  are  .06tt  seconds  apart.  Figure  2-2-10  shows  the  noise 

C 

corrupted  output,  y(t)  + n(t),  for  /N  - 10 dB. 

The  system  is  first  identified  for  n(t)  = 0.  The  operator  generating 
the  set  of  basis  functions  is  chosen  to  be  the  advance  operator  discussed 
earlier. 


Case  1 - Zero  Additive  Noise  (n(t)  = 0) 


The  input  and  output  basis  functions  generated  for  an  advance  of  10 
time  samples  (i.e.,  .6ir  seconds)  are  plotted  in  Figures  2-4-11  and  2-4-12 
respectively. 


When  a least-squares  approximation  scheme  is  used,  the  homogeneous 


set  of  equations  is  represented  in  matrix  form  by 


(2-4-76) 


0 


000  ?-30  4.00  6.00  a. 00  10.00  12.00  14.00  16.00  1 8.00  20.00 

Fig.  2-4-9.  Pseudo-random  uniformly  distributed  noise  sample,  n(t). 


Fig.  2-4-11.  Basis  set  for  the  input. 


.00 
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The  first  n+1  columns  of  [M]  represent  the  y(.)  basis  functions,  and  the  last 
nrt-1  columns  of  [M]  represent  the  x(.)  basis  functions. 

The  poles  and  zeros  resulting  from  the  identification  are  presented  in 
Table  2-4-7  along  with  the  exact  poles  and  zeros.  Clearly  the  discrepancy 
between  the  exact  and  identified  values  is  due  to  the  numerical  convolution  that 
generated  y(t). 


Table  2-4-7 


Poles 

Zeros 

Exact 

Identified 

Exact 

Identified 

real 

imaginary 

real 

imaginary 

real 

imaginary 

real 

imaginary 

-.3183 

.5 

-.3173 

.5027 

-.4411 

0 

-.4305 

0 

-.3183 

-.5 

-.3173 

-.5027 

-.3478 

1.0086 

-.3816 

1.013 

-.1592 

1 

-.1602 

1.001 

-.3478 

-1.0086 

-.3816 

-1.013 

-.1592 

-1 

-.1602 

-1.001 

— 

— 

Case  2 - Nonzero  Additive  Noise 

When  the  output  is  corrupted  by  additive  noise  the  set  of  basis  functions 
generated  by  the  advance  operator  of  Case  1 is  given  by 

{y(t  + .6irk)  + n(t  + .6irk)};  k- 0,1,2,3,*<.  (2-4-77) 

The  set  of  basis  functions  of  (2-4-77)  is  plotted  in  Figure  2-4-13.  Since  the 
output  y(t)  is  a decaying  signal  and  the  noise  n(t)  is  uniformly  distributed, 
the  is  time  dependent.  The  S/N  associated  with  each  basis  function  is 
presented  in  Figure  2-4-13. 


j 
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It  is  clear  that  in  practice  the  output  y(t)+n(t)  cannot  be  separated 
into  y(t)  and  n(t) • But,  for  the  given  synthesized  example  we  will  define 
their  contribution  through  separate  matrices  in  order  to  gain  insight  into  the 
noise  effect.  Define  the  noise  matrix  [N]  such  that  its  first  n+1  columns 
are  the  advanced  versions  of  the  noise  sample  n(t),  and  the  last  m+1  columns 
are  zero  entries  (i.e.,  no  noise  associated  with  the  input).  In  particular 


[N] 


I 


I 


0 


I 


I 


n(t.)  I n(t . + .6  ) | • • .ln(t.+2.4  ) I 

1 . l 1 .1  i 


I 


n+1 


-1°: 


i o 


(2-4-78) 


acr 


For  the  noisy  case,  (2-4-76)  becomes 


( [M]T  + [N]T  } { [M]  + [N]  } 


D„ 


A 

N„ 


L N_ 


0 

0 


0 J 


(2-4-79) 


From  (2-4-79)  it  is  found  out  that  the  noise  contribution  is  manifested  by 


[N^]  defined  by 


[N  1 - ([M]T[N]  + [N]T[M]  + [N]T[N]}. 


(2-4-80) 


■t-Jarntl- 


1 
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Carrying  out  the  matrix  operations  in  (2-4-80)  results  in 

n+1  nri-1 


■ 


n+l 


nrt-1 


2<n(.),y(.)>  + <n(.)»n(.)>  j <n(.),x(.)> 

• (2-4-81) 

<n( . ) , x(.)>  , 0 


The  inner  product  matrix  for 


Case  1 (without  noise)  is  given  by 


[M]T[M] 


ISO , 3 
21.08 
*70.95 
1.418 
24.43 
18.6 
"6.466 
"2.208 
4.992 


21.08 

99.16 

"12.12 

"52.55 

21.19 

24.78 

6.454 

"7.374 

0.734 


"70.95 

1.418 

26.43 

18.6 

"6.466 

"2.280 

4.992 

“12.12 

"52.55 

21.19 

24.78 

6.454 

"7.374 

0.734 

43.87 

"1.225 

"19.04 

"9.61 

4.498 

1.59 

"3.251 

"1.225 

32.77 

"7.307 

"13.42 

"4.95 

4.001 

0.4380 

"19.04 

"7.307 

12.62 

6.506 

"1.087 

"1.902 

1.542 

"9.61 

"13.42 

6.506 

8.553 

1.254 

"1.963 

0.5065 

4.498 

"4.95 

"1.087 

1.254 

1 + 248 

"0.4067 

"0.4057 

1.59 

4.001 

"1.902 

"1.963 

"0.4067 

0.5029 

“0.09747 

"3.251 

0.4300 

1.542 

0.5065 

"0.4057 

“0.09747 

0.2655 

(2-4-82) 


For  * 10  dB,  the  contribution  due  to  the  noise  is 


(Ntl  - 


13.72 

“8.070 

"5.345 

4.132 

3.719 

"0,3761 

"0.4894 

0.1647 

0.2172 


"8.070 
13.13 
"2.664 
"5 • 029 
4.806 
0.1295 
0.368 
"0.2664 
"0.09355 


“5.345 

"2.664 

14.19 

"3.339 

"2.596 

“0.5919 

0.4912 

0.2192 

"0.3601 


4.132 

“5.829 

"3.339 

13.56 

“4.261 

"1.015 

"0.6522 

0.4182 

0.1588 


3.719 

4.806 

"2.596 

"4.261 

0.553 

1.516 

0.2575 

"0.4491 

0.07427 


"0.3761 
0. 1295 
"0.5919 
"1.015 
1.516 
0 
0 
0 
0 


"0.4894 

0.368 

0.4912 

"0.6522 

0.2575 

0 

0 

0 

0 


0.1647 

"0.2664 

0.2192 

0.4182 

"0.4491 

0 

0 

0 

0 


0.2172 

"0.09355 

"0.3601 

0.1588 

0.07427 

0 

0 

0 

0 


(2-4-83) 


Then,  solving  the  equation 


([M]T[M]  ♦ [Nj} 


(2-4-84) 


1 
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yields  the  hatted  coefficients  for  the  characteristic  equation,  whose  roots, 
after  the  proper  transformation  result  in  the  identified  poles  given  by 


.0464  + j 0.689 
.165  + j 1.24  . 


(2-4-85) 


From  Figure  2-4-13,  it  is  clear  that  there  exists  a trade-off  between 

the  amount  of  advance  used  and  the  given  S/N  . in  this  respect,  the  present 

example  does  not  represent  any  effort  to  optimize  the  advance  used. 

To  improve  the  solution  obtained  in  (2-4-85) , one  should  1)  increase 

the  number  of  samples  used  over  the  given  record  (the  basis  functions  of  the 

example  utilized  just  61  samples) , 2)  properly  choose  a usable  portion  of  the 

record,  and  3)  choose  a 'good'  value  for  the  amount  of  advance.  In  addition, 

the  approximation  scheme  can  be  improved  by  taking  advantage  of  the  noise 

properties.  For  instance,  in  a least-squares  scheme  we  premultiply  the 

equation  by  the  transpose  of  the  rectangular  matrix  representing  the  noisy 

basis  functions  (see(2-4-79) ) . Define  the  matrix  operation  0 as  first 

reversing  the  order  of  the  matrix  rows,  and  then  taking  the  transpose. 

Utilizing  this  operator  in  (2-4-79)  instead  of  the  regular  transpose  might 

reduce  the  effect  of  noise.  This  operation  takes  advantage  of  the  small 

correlation  between  noise  samples  separated  by  large  time  intervals.  In 

T 

particular,  one  contributer  to  [N  ] in  (2-4-80)  is  [N]  [N] . Concentrating 

L 

T 

on  the  nonzero  submatrix  of  [N]  [N]  of  Case  2,  we  have 


[N]T[N] 


3.266 

'4.13 

"4.924 

3.778 

2.959 


1.566 

3.327 

'1.989 

'5.078 

6.018 


'1.832 

1.45 

4.062 

'1.714 

'1.159 


2.401 

'2.36 

'0.06135 

3.675 

'3.034 


0.6318 

0.5546 

'2.535 

'0.3168 

1.264 


(2-4-86) 


The  same  matrix  with  a 0 operation  results  in 


1.824  ' 

“1.309 
0.02368 
“0.2852 
“0,1529 


(2-4-87) 


Notice  that  21  out  of  25  entries  of  (2-4-87)  are  smaller  in  magnitude  than 
the  corresponding  entries  in  (2-4-86). 


[N]G[N]  = 


.2258 
0.1271 
0.1849 
“0.1706 
. 1 . 824 


0.1271 

0.4922 

"0.8318 

0.9663 

"1.309 


0.1849 
"0.8318 
“0 . 2702 
“2.058 
0 . 02368 


"0.1706 
0.9663 
"2.058 
0.7081 
"0 . 2852 


This  reduces  the  effect  of  the  noise. 
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USEFUL  SUPPLEMENTAL  TOOLS 


Two  supplemental  tools  are  discussed  in  this  chapter.  The  first 
is  the  concept  of  the  "growing  tree."  The  "growing  tree"  describes 
the  movement  of  poles  in  the  S-plane  as  the  order  of  the  approximation 
is  increased.  The  tree  behaves  in  a distinct,  but  highly  structured 
manner  for  each  signal.  Consequently,  the  tree  is  useful  for  extra- 
polating the  poles  of  a lower-order  approximation  to  those  of  higher 
order.  Also,  the  tree  is  useful  in  separating  signal  poles  from 
extraneous  poles  in  a noisy  problem. 

The  second  topic  deals  with  the  problem  of  obtaining  a "good" 
first  guess  for  iterative  schemes  which  utilize  highly  nonlinear 
equations. 

3.1  The  Growing  Tree 

Assume  an  rth-order  approximation  is  to  be  used  to  approximate  a 
function  which  is  composed  of  a sum  of  n exponentials.  If  the  order 
of  approximation  is  increased  from  r=l  to  r=n,  a trajectory  of  the  pole 
locations  can  be  drawn  in  the  S-plane  (see  Fig.  2-4-6).  Clearly,  for 
r>n,  the  set  of  equations  for  determining  the  coefficients  becomes 
singular.  Hence,  once  the  exact  solution  has  been  reached,  corre- 
sponding to  r-n,  the  "growth"  of  the  trajectory  is  stopped. 

Now  assume  that  a nonexponential  function  is  to  be  approximated 
by  a sum  of  exponentials.  From  Chapter  1,  it  is  known  that  an  infinite 
set  of  exponential  basis  functions  spans  the  space  L^fo.00]*  Hence,  if 


3-2 


1 


the  nonexponential  function  belongs  to  L^tO,00],  the  function  can  be 
approximated  arbitrarily  closely  by  the  exponential  basis.  Since  it 
is  impractical  to  use  an  infinite  set  as  a representation,  an  approx- 
imation of  order  r is  utilized.  A trajectory  of  the  pole  locations  is 
structured  in  the  S-plane  as  r is  increased  successively  from  unity. 

This  "growing"  trajectory  is  called  a "growing"  tree.  Clearly,  the 
tree  of  a nonexponential  function  continues  to  grow  with  ever-increasing 
r. 

When  discussing  the  behavior  of  the  growing  tree,  it  is  convenient 
to  divide  the  nonexponential  functions  into  two  subgroups.  A natural 
division  consists  of  1)  functions  which  have  a faster  than  exponential 
decay,  and  2)  functions  which  have  a slower  than  exponential  decay.  The 
first  group  is  considered  to  be  composed  of  time-limited  functions 
(in  some  wide  sense)  while  the  second  group  is  considered  to  be  com- 
posed of  nontime-limited  functions. 

On  the  basis  of  many  worked  examples,  it  is  noticed  that  the 
growth  of  a tree  belonging  to  a time-limited  function  expands  in  the 
general  direction  of  the  negative  real  axis  of  the  S-plane.  Also,  if 
r^r^,  the  poles  representing  r^  are  positioned  on  a contour  of 
greater  extent  than  the  contour  of  poles  representing  r^.  The  two 
contours  do  not  intersect.  For  real  time  signals,  the  tree  is  symmet- 
rical about  the  negative  real  axis.  As  r is  increased,  the  tree 
spreads  out  in  the  positive  and  negative  directions  of  the  imaginary 
axis  as  it  grows  along  the  negative  real  axis. 

For  nontime-limited  functions,  the  base  of  the  tree  is  localized 


in  a narrow  region  of  the  negative  real  axis  with  the  branches  extending 
towards  the  imaginary  axis.  Since  the  functions  being  approximated 


are  assumed  to  decay  with  increasing  time,  the  poles  do  not  cross  over 
into  the  right-half  plane  as  r is  increased  in  value.  Instead  the 
branches  stop  growing  just  short  of  the  imaginary  axis.  For  real 
signals  the  tree  is  again  symmetrical  about  the  negative  real  axis. 
However,  because  nontime-limited  functions  are  bandlimited , the  tree 
branches  are  confined  within  a finite  region  of  the  imaginary  axis. 

For  r2>r^>  the  contour  formed  for  r^  grows  inside  that  for  r^  by 
interleaving  with  all  of  the  contours  corresponding  to  r<r  . Hence, 
the  branches  of  the  tree  tend  to  bend  inwards  (as  though  they  were 
heavily  laden  with  fruit). 

Some  examples  are  now  presented  to  illustrate  the  behavior  ot  the 
two  types  of  growing  trees. 

Example  3-1-1 

Consider  a strictly  time-limited  sine  pulse,  denoted  by  S1NP,  which 
is  given  by 

(sinirt  ; O^t^l 

(3-1-1) 

0 ; elsewhere. 

The  SINP  tree  with  r = 1,  2,  . . . , 8 is  presented  in  Fig.  3-1-1.  The 
SINP  function  and  its  approximations  for  odd  and  even  values  of  r, 
respectively,  are  presented  in  Figures  3-1-5  and  3-1-6.  The  integrated 
squared  error,  as  a function  of  r,  is  plotted  in  Fig.  3-1-11. 
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Example  3-1-2 

Consider  a strictly  time-limited  square  pulse,  denoted  by  SQP,  which  is 
given  by 


(1;  0<t<l 
0;  elsewhere  . 


('i-  1-2) 


The  SQP  tree  with  r = 1,  2,  . . . , 8 is  presented  in  Fig.  3-1-2.  The 
SQP  function  and  its  approximations  for  odd  and  even  values  of  r, 
respectively,  are  presented  in  Figures  3-1-7  and  3-1-8.  The  integrated 
squared  error,  as  a function  of  r,  is  plotted  in  Fig.  3-1-11. 


Example  3-1-3 

Consider  a wide  sense  time-limited  gaussian  pulse,  denoted  by  GP,  which 
is  given  by 


GP 


r -[4(t-i/2)]' 

e 


t£0 
t<0  . 


(vt-*) 


The  GP  tree  with  r = 1,  2,  . . .,  8 is  presented  in  Fig.  3-1-3.  The 
GP  function  and  its  approximations  for  odd  and  even  values  of  r, 
respectively,  are  presented  in  Figures  3-1-9  and  3-1-10.  The  integrated 
squared  error,  as  a function  of  r,  is  plotted  in  Fig.  3-1-11. 


In  carrying  out  the  approximations  for  the  above  three  examples, 
the  reverse  time  integral  operator  (see  (2-4-43)  ) was  used  along  with 
the  ISE  approximation  scheme  (see  (A-l-11)  ).  For  the  final  example 
the  differential  operator  is  utilized  along  with  the  ISE  approximation 


scheme. 
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Example  3-1-4 

Consider  the  nontime-limited  function,  denoted  by  SINC,  which  is  given  by 


SINC 


/ sin_t 

”t"  ; t*0 

k 0 ; t<0  . 


TEXT  CONTINUED  FOR  THIS  EXAMPLE  ON  PAGE  3-15. 


'O.CO  2.50  5 OO  7.  btl  10.00  12.50  15.00  17.50  20.00  22.50  25  00 

uio-i) 

Fig.  3-1-6.  Even  order  approximations  of  SINP. 


I 


'c.ro  'so  5.uc  7'.  o,i>  ;p.nc  :o.so  ib.ro  iV.  io  FcrTcc  o’o.so  ob.oo 
(XIO-iJ 

Fig.  3-1-8.  Even  order  Approximations  of  SQP. 


? -t 1 1 1 1 — 1 1 I 1 1— 

0.00  2.50  5.  CO  7.  SO  10.  CO  12.  SO  15.00  17.  SO  20.00  22.50  25.00 


U10-!) 

fig-  3-1-9.  Odd  order  approximations  of  GP. 
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Although  this  function  exists  for  all  positive  time,  only  that  portion 
over  the  interval  [0,  4nt]  was  approximated  (i.e.,  the  integrations  in 
the  inner  products  extended  from  0 to  41t) . The  SINC  tree  with 
r = 1,  2,  . . . , 8 is  presented  in  Fig.  3-1-4.  The  SINC  function  and 
its  approximations  for  odd  and  even  values  of  r,  respectively,  are 
presented  in  Figures  3-1-12  and  3-1-13.  The  integrated  squared  error, 
as  a function  of  r,  is  plotted  in  Fig.  3-1-14. 

Notice  that  the  "time-limited"  trees  in  Figures  3-1-1  through 
3-1-3  grow  to  the  left  whereas  the  "nontime-limited"  tree  in  Fig.  3-1-4 
grows  to  the  right.  The  trees  demonstrate  that  a priority  exists  for 
the  manner  in  which  a waveform  is  approximated  as  the  number  of  poles 
is  increased.  Given  a small  value  of  r,  the  poles  are  chosen  so  as 
to  give  a "best"  fit  to  that  portion  of  the  waveform  where  most  of  the 
energy  is  concentrated.  Additional  poles  are  then  used  to  refine  the 
approximation  to  the  remaining  portion  of  the  signal.  For  a time- 
limited  signal  the  additional  poles  add  terms  to  the  approximation 
which  have  successively  larger  rates  of  decay  in  an  effort  to  account 
for  the  faster-than-exponential  decay  of  the  signal  at  later  instants 
of  time.  Conversely,  for  a nontime-iimited  signal  the  additional 
poles  add  terms  to  the  approximation  which  have  successively  slower 
rates  of  decay  in  an  effort  to  account  for  the  slower-than-exponential 
decay  of  the  signal  at  later  instants  of  time. 

For  smooth  time  functions  the  tree  branches  consist  of  relatively 


smooth  contours  which  are  U-shaped  about  the  real  axis.  The  faster  the 
rate  of  growth  along  the  real  axis,  the  faster  is  the  decrease  in  the 


Fig-  3-1-12.  Odd  order  approximations  of  SINC. 
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Fig.  3-1-14 


The  ISE  as  a function  of  r for  SINC 
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ISE.  On  the  other  hand,  piecewise-continuous  functions  have  tree 
branches  which  tend  to  wiggle  to  the  right  and  left  during  the  growth 
of  the  tree.  Since  these  branches  do  not  progress  very  rapidly  along 
the  real  axis,  their  ISE  decreases  more  slowly.  This  is  illustrated 
by  the  ISE curves  of  Figures  3-1-11  and  3-1-14.  This  behavior  is 
manifested  in  the  time  waveform  approximations  as  a Gibb's  phenomenon 
(see  Figures  3-1-5  through  3-1-8). 

The  growth  of  the  trees  is  seen  to  be  highly  structured.  Therefore, 

given  a tree  up  to  order  r,  it  is  possible  to  obtain  a reasonably  good 

guess  of  the  (r  + l)1"*1  trajectory  by  simply  extrapolating  each  tree 

branch.  In  Fig.  3-1-2  for  the  SQP  tree,  the  extrapolated  pole  positions 

for  r = 9 are  shown  along  with  the  actual  approximation. 

t h 

It  is  interesting  to  observe  that  the  r order  trajectory  for 
both  time-limited  and  nontime-limited  functions  is  U-shaped  to  the 
right  about  the  real  axis.  The  only  difference  is  that,  for  time- 
limited  functions,  the  dominant  poles  are  at  the  tips  of  the  U while, 
f> r nontime-limited  functions,  the  dominant  poles  are  at  the  base  of 
the  U.  It  is  also  interesting  to  note  that,  although  the  trees  in 
Figures  3-1-1,  3-1-2,  and  3-1-3  are  for  three  different  pulses  having 
approximately  the  same  pulse  width,  the  trees  span  roughly  the  same 
region  of  the  S-plane. 

The  growing  tree  is  a useful  tool  for  discriminating  against  noise. 

For  a signal  composed  of  n exponentials  the  tree  stops  growing  at  r = n. 
When  noise  is  present,  the  tree  continues  to  grow  beyond  r = n.  However, 
the  actual  poles  tend  to  be  perturbed  only  slightly  while  the  additional 
poles,  due  to  the  noise,  move  significantly.  This  behavior  can  be  used 
to  determine  both  the  system  order  and  system  poles  under  noisy  conditions. 
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3 . 2 Iterative  Optimization  Techniques 

The  approximation/identif ication  technique  discussed  previously  are 
linear  in  the  sense  that  a linear  set  of  equations  are  solved  for  the 
numerator  and  denominator  coefficients  of  the  transfer  function.  These 
equations  are  generated  by  selecting  the  coefficients  to  minimize  an 
error  criterion  such  as  least  squares,  integrated  squared-error,  etc 
The  system  poles  and  zeros  are  nonlinearly  related  to  the  denominator  and 
numerator  coefficients,  respectively,  by  means  of  polynomial  equations. 
Some  approximation/identif ication  techniques  solve  directly  for  the  poles 
and  residues.  The  equations  are  now  generated  by  selecting  the  poles  and 
residues  to  minimize  some  suitable  error  criterion.  These  equations  are 
highly  nonlinear  and  are  usually  solved  through  some  iterative  scheme. 
Different  iterative  methods,  mainly  designed  for  the  synthesis  problem, 
are  presented  in  [19],  [20],  [21],  [22],  [27],  [28].  A good  initial 
guess  is  highly  desirable  in  order  to  minimize  the  number  of  iterations. 
In  some  cases,  a good  initial  guess  is  required  if  the  correct  solution 
is  to  be  achieved  due  to  the  presence  of  multiple  local  minima  (i.e.,  the 
iteration  should  begin  near  the  correct  local  minimum) . It  should  be 
noted  that  minimizing  an  error  criterion  with  respect  to  the  transfer 
function  coefficients  is  not  equivalent  to  minimizing  an  error  criterion 
with  respect  to  the  poles  and  residues.  However,  both  approaches  tend 
to  be  equivalent  when  the  errors  are  small. 

Although  any  one  of  the  identification  techniques  discussed  in 
chapter  2 could  be  used  to  provide  the  initial  guess,  not  all  of  the 


L 


operators  and  approximation  techniques  perform  the  same  for  a given  pro- 
blem. Hence,  the  identification  technique  chosen  for  a specific  problem 
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should  be  selected  with  care.  This  is  illustrated  by  the  example  dis- 
cussed next. 

Example  3-2-1 

This  example  is  one  of  the  problems  given  at  the  RADC  Spectral  Es- 
timation Workshop  held  May  24-26,  1978.  It  is  interesting  because  it 
involves  an  FM  signal  which  can  be  interpreted  in  terms  of  a slowly 
varying  linear  single-tuned  circuit. 

This  problem  concerns  identifying  emitters  by  their  signal  character- 
istics . 

Problem:  Determine  the  instantaneous  frequency  of  the  largest  am- 
plitude sinusoid  within  the  portion  of  the  spectrum  of  interest.  The 
instantaneous  frequency  is  defined  as 

2tt  dt 

where  <{>  = argument  of  the  sinusoid. 

It  is  required  to  determine  the  instantaneous  frequency  every  0.04  sec. 

from  t,  through  t,..  where 
1 16 

t^  » 0.04  sec. 
t2  * 0.08  sec. 

t,,  = 0.64  sec. 

16 

The  data  model  is 

P 

x(t)  cos  (0^)  + n(t)  (where  P £ 5) 

The  summation  consists  of  a single  large-amplitude  sinusoid  (the  "signal") 
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plus  a number  of  smaller-amplitude  sinusoids  (the  "interference") . There 
are  "P"  sinusoids  in  all. 


The  "noise"  waveform  is  n(t).  The  resulting  data  waveform,  x(t), 
has  been  sampled  at  an  800  samples/sec.  rate,  and  the  total  duration  of 
the  data  waveform  is  0.64  sec.  The  sampled  noise  waveform  is  actually 
a sequence  of  pseudo-uniformly  distributed  random  numbers. 

The  spectrum  of  interest  is  from  130  to  270  Hertz.  The  sampling 
rate  was  chosen  to  be  at  least  four  (4)  times  greater  than  the  instan- 
taneous frequency  of  any  sinusoid  in  the  summation  above. 

Some  statistics  concerning  the  data  are: 

1.  The  rate  of  frequency  drift  of  any  sinusoid  is  not  greater  than 
190  Hertz/sec. 

2.  An  estimate  of  the  signal-to-nolse  ratio  is: 


where  N = 512  = the  total  number  of  data  samples, 

denotes  the  sampled  value  of  the  largest  amplitude  sinusoid, 
without  noise,  at  the  time  point  t^ 

n^  denotes  the  sampled  value,  at  time  point  t^,  of  the  noise  waveform, 
n(t)  . 


3.  An  estimate  of  the  signal-to-noise-plus-interference  ratio  is 
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N 2 

l Si 

1=1  <\» 

10  log1Q 17  dB 

N P-1 

(ni  + l 

1=1  K-l 

where  and  n^  have  already  been  defined  in  (2)  and  where 
P-1 

£ 5C.  denotes  the  sum  of  the  sampled  values,  at  time  point  t , of  all 

K=1  Ki  1 

of  the  P-1  weaker  sinusoids,  the  "interference". 

The  first  thing  to  be  noted  in  this  problem  is  the  low  sampling 
rate  which  may  produce  as  few  as  four  samples  per  cycle.  In  addition, 
the  frequency  drift  causes  the  problem  to  be  time  variant.  If  the  pro- 
blem is  to  be  interpreted  in  terms  of  a sequence  of  time-invariant  sys- 
tems, only  a limited  number  of  samples  can  be  used  to  estimate  the  fre- 
quency drift  rate  of  190  Hz/sec.,  the  maximum  frequency  drift  in  .04 
seconds  in  7.6  Hz.  Since  this  is  less  than  6%  of  the  lowest  frequency 
in  the  spectrum  of  interest  (i.e.,  130  Hz),  the  instantaneous  frequency 
is  assumed  to  be  approximately  constant  over  this  interval.  Hence,  32 
samples,  spaced  about  the  data  point,  were  used  to  determine  the  instan- 
taneous frequency  at  each  instant  t^.  The  slow  sampling  rate  dictates 
the  use  of  a discrete  time  operator  for  the  identification  procedure 
(i.e.,  it  is  difficult  to  visualize  the  time  waveform  from  the  data). 

In  addition,  the  small  number  of  samples  which  can  be  used  at  a time, 
discourages  the  use  of  a delay  operator.  Consequently,  it  was  decided 
to  use  an  overdetermined  Prony  approach.  Because  no  exponential  decays 
are  involved  in  this  problem,  ill  conditioning  due  to  extremely  small 
sample  values  is  not  a concern. 
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Since  complex  poles  in  a real  system  appear  in  conjugate  pairs,  the 
overdetermined  Prony  method  was  used  to  solve  for  two  poles.  (The  inter- 
ference and  noise  were  assumed  to  be  negligible  since  the  signal-to-noise- 
plus-interference  ratio  was  given  to  be  approximately  17  dB) . These 
poles  provided  the  initial  guess  for  a least  squares  direct  search 
iterative  scheme  which  varied  the  poles  and  residues  in  a systematic 
manner  in  order  to  minimize  the  error  between  the  approximation  and  the 
sample  values  in  the  .04  second  interval  being  considered. 

The  exact  solution  to  the  synthesized  problem,  as  provided  by  RADC, 
is  presented  in  Appendix  B.  The  initial  guess  by  means  of  the  Prony 
technique  and  the  solution  obtained  via  the  direct  search  iteration  are 
shown  in  Figure  3-2-1  along  with  the  exact  solution  (solid  line).  Fig. 

3-2-2  is  a plot  of  the  instantaneous  error  for  both  the  initial  guess 
and  the  direct  search.  The  initial  guesses  are  seen  to  be  quite  good. 

However,  the  iteration  scheme  does  result  in  an  improvement  at  13  of  the 
16  points.  This  solution  was  the  best  of  all  those  submitted  at  the 
workshop  and  was  considerably  better  than  those  obtained  by  conventional 
spectral  estimation  techniques  such  as  Burg's  maximum  entropy,  FFT,  and 
ARMA  (autoregressive  moving  average  filter)  approaches.  A technique, 
based  upon  zero-crossing  FM  demodulation  [33],  provided  a solution  which 
was  practically  as  good.  Results  for  this  approach  are  also  plotted  on 
Fig.  3-2-1  and  compared  in  Fig.  3-2-3  where  the  zero-crossing  method  is  seen  to 
provide  closer answe . s in  6 of  the  16  points.  The  sum  of  squared  errors 
for  the  iterative  approach  is  2.34  whereas  that  for  the  zero-crossing 
approach  is  4.19. 
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Fig.  3 2-2.  Instantanuous  errors  for  the  Prony/direct  search 
method . 
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Chapter  4 


IDENTIFICATION  OF  WEAKLY  NONLINEAR  SYSTEMS 


4.1  Introduction 

The  Volterra  series  representation  for  a weakly  nonlinear  system  was 
introduced  in  Chapter  1.  The  series  is  most  useful  in  applications  where 
a finite  number  of  terms  N is  adequate  for  modeling  the  system  response 
(see  (1-1-6) ) . It  has  been  shown  [34]  that  the  nth-order  impulse  response 
;i > • • • > t ) is  a sum  of  exponentials  when  the  linear  impulse  response 
h^(t)  is  a sum  of  exponentials.  Also,  the  natural  frequencies  of  the 
higher  order  transfer  functions  can  be  determined  a priori  from  the  natural 
frequencies  of  h^(t). 

In  particular,  if  the  linear  impulse  response  is  given  by 


hx(t) 


K , 

.1  \ * S 

1=1  ; t > 0 


0 ; t < 0, 


(4-1-1) 


then  the  second-order  impulse  response  is 

/ M K 
k 


h2(tl,t2) 


ll  \k  e*i 

ij-1  k2=l  *1*2 


C1  + \2t2 


; c2  - ci 


(4-1-2) 


re  \ C2  + ak  C1 

l l \ k e 1 2 • t < t 

.k^l  k2=l  12  ’ C2  - 1 
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The  terms  for  the  third-order  impulse  response  are  of  the  same  form  as  in 

(4-1-2).  However,  now  a triple  summation  is  involved  over  the  indices 

k^ , k9,  and  k^  and  the  exponents  of  each  term  are  of  the  form 

(a  • t + a - t + a • t ) where  i,  j,  l = 1,2,  or  3,  i / j ^ SL,  and  there 
i 1 Kj  Z *4  3 

are  six  terms,  one  for  each  condition  t.  >_  t.  >_  t.  (See  (7-5)  of  [34]). 

i J 

The  relation  between  the  natural  frequencies  of  (4-1-2)  and  those  of 
(4-1-1)  is  given  by 


(4-1-3) 


In  addition,  the  parameter  M in  (4-1-2)  is 

M = K2  + 1 . 


(4-1-4) 


Hence,  given  h (t),  the  only  unknowns  in  (4-1-2)  are  the  residues  A , . 

1 kl*2 

These  are  determined  using  an  appropriate  sum  of  exponentials  for  the 

input  to  the  weakly  nonlinear  system.  The  response  is  a sum  of  exponentials. 

This  response  is  assumed  to  be  the  output  of  an  equivalent  linear  system. 

The  unknowns  A , are  unique-lv  related  to  the  residues  of  the  equivalent 
k12 

linear  impulse  response.  Hence,  by  applying  one  of  the  approximation 

schemes  to  the  equivalent  impulse  response,  a set  of  linear  equations  can 

be  written  from  which  the  coefficients  A^  , can  be  determined. 

*12 

The  procedure  begins  by  first  determining  h^(t).  For  this  purpose, 
the  input  amplitude  is  maintained  small  enough  such  that  nonlinear  effects 
are  negligible.  The  next  step  is  to  determine  ^2^1’ t2^’  7,16  input 


4-3 


amplitude  is  now  chosen  such  that  second-order  effects  are  noticeable 
but  third-and  higher-order  effects  are  still  negligible.  The  response  due 
to  the  second-order  portion  of  the  system  can  be  isolated  by  subtracting 
out  a suitably  weighted  replica  of  the  linear  response.  The  procedure 
continues  for  h^t^.t^.t^)  and  so  on.  To  determine  the  appropriate  signal 
amplitude  at  each  step,  a single  sinusoid  is  used  as  the  input  and  the 
harmonic  content  of  the  output  is  checked. 

4 . 2 Identification  of  a Transistor  Amplifier 

The  identification  approach  discussed  above  was  applied  in  [34]  to 
identification  of  a transistor  amplifier.  The  example  is  repeated  here 
where  the  results  of  this  dissertation  are  utilized  to  give  additional 
insight  into  the  problem.  The  amplifier,  with  the  transistor  replaced  by 
its  nonlinear  incremental  equivalent  circuit,  is  shown  in  Fig.  4-2-1. 


son 


Fig.  4-2-1.  Transistor  Amplifier 
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t 

The  circuit  and  model  parameters  are  given  in  [34]  and  are  not  repeated 

I 

here.  From  [34]  it  is  known  that  the  system  behaves  linearly  for  input 

-3  -2 

amplitudes  less  than  10  volts.  With  an  input  amplitude  of  10  volts, 

a second-order  response  is  noted  while  higher-order  responses  remain 

negligible. 

In  [34]  the  first-order  impulse  response  was  identified  using  as  an 

input 


v.  (t) 
in 


.001 


-107t 


t > 0 


t < 0 


(4-2-1) 


Identification  of  the  second-order  impulse  response  utilized  the  input 


r 


; t •>  0 

; t < 0 


(4-2-2) 


The  second-order  portion  of  the  response  was  separated  from  the  total 
response  (i.e.,  the  sum  of  first-and  second-order)  by  subtracting  from  the 
total  response  10  times  the  linear  response. 

The  identification  procedure  used  in  [34]  is  Jain's  pencil  of  functions 

technique  [31].  The  operator  in  this  approach  is  the  simple  integral 

t 

operator,  / (•)  dT.  In  order  to  obtain  meaningful  results  with  this 

o 

operator,  it  was  necessary  to  use  2400  samples  closely  spaced  at  intervals 
-9 

of  10  seconds.  The  system  order  was  established  to  be  n = 2.  (This 
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value  is  correct,  as  is  obvious  by  inspection  of  Fig.  4-2-1).  The  poles 
and  residues  found  in  [34]  for  h^(t)  are  given  in  Table  4-2-1. 


Table  4-2-1 


residue 


-.011551(2-10  tt) 


.2807(10  ) 


-10.6169(2-10  tt) 


273.684(10  ) 


As  pointed  out  in  case  5 of  example  2-4-1,  the  simple  integral  operator 
is  a poor  choice  because  it  results  in  a highly  linearly  dependent  set  of 


basis  functions.  Consequently,  the  example  was  repeated  using  the  advance 

(negative  delay)  operator  of  (2-4-39) . Only  151  samples  spaced  at  intervals 
_8 

of  10  seconds  were  used  (i.e.,  many  fewer  samples  and  a slightly  shorter 
data  record  than  in  [34]).  The  advance  T was  set  equal  to  10  2 seconds 
(i.e.,  every  10  samples).  The  poles  and  residues  determined  in  this 
manner  are  presented  in  Table  4-2-2. 


Table  4-2-2 


residue 


= -.011454(2-10  tt) 
X2  = -10.8071  (2- 106tt) 


.27502(10  ) 
279.107 (106) 


Figure  4-2-2  presents  the  true  output  and  the  two  responses  generated  by 
the  poles  and  residues  of  Tables  4-2-1  and  4-2-2. 

From  (4-1-2)  notice  that  h2(t^,t2),  evaluated  for  = t2  = t,  becomes 


M K 


h.(t,t)  * h equlv.  (t)  - J ^ ^ 1c  G * 


(4-2-5) 


From  (4-1-4)  and  table  4-2-1,  it  follows  that 


M=K2+1=22+1=5  . 


IX10-7) 

Fig.  4-2-2.  True  output  and  responses  generated  by  the  identified  poles 
and  residues. 


4-7 


Hence,  the  natural  frequencies  of  (4-2-3)  are  tabulated  in  Table  4-2-3. 

Table  q-2-i 


# 

M 

K 

Natural  frequency 
of  (4-2-3) 

in  terms  of  \2 

Camments 

1 

1 

1 

31  + ai 

2\ 

2 

2 

1 

a2  + al 

X2  + \ 

3 

3 

1 

a3  + ai 

2A1  - A2 

Not  allowed 

4 

4 

1 

^l 

A2 

5 

5 

1 

a5  +al 

A! 

6 

1 

2 

al  + a2 

A1  + A2 

Same  as  2 

7 

2 

2 

a2  +a? 

“*r 

8 

3 

2 

1 2 

A! 

Same  as  5 

9 

4 

2 

— 

S4  +S? 

2A2  - \ 

Not  allowed 

10 

5 

2 

35  + 3 2 

A 

l2 

Same  as  4 

Therefore,  the  allowable  natural  frequencies  of  (4-2-3)  are  given  by 

*2»  2A1,  2A2,  A1  + A2}  . (4-2-5) 

A procedure  for  determining  the  coefficients  A^  ^ is  outlined  in  [34], 


It  should  be  mentioned  that  the  set  of  equations  for  A.  becomes 

lk2 

highly  linearly  dependent  when  the  dynamic  range  in  the  natural  frequencies 


of  the  linear  response  is  large.  For  example,  if  A2  >>  A^,  then 


A2  + Aj  r A2>  The  natural  frequencies  of  (4-2-5)  then  give  rise  to  an 
almost  linearly  dependent  set  of  basis  functions.  In  our  example. 


10' 


(4-2-6) 
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Consequently,  the  equations  in  [34]  for  A,  , were  extremely  difficult 

kl  2 

to  solve.  A remedy  for  this  situation  is  to  consider  A^  + A^  to  be  the 
same  natural  frequency  as  A^.  As  a result,  the  set  of  (4-2-5)  is  reduced 


i\,  X2,  2Ai,  2A2) 


(4-2-7) 


Removal  of  the  natural  frequency  (A^  + A2)  eliminates  the  ill-conditioned 
nature  of  the  equations  without  severely  affecting  the  identification. 


I 


Chapter  5 


SUMMARY 

5.1  Discussion  of  Results  and  Conclusions 

A general  parametric  model  for  the  black  box  identification  of 
linear  time  invariant  systems,  in  terms  of  rational  transfer  functions, 
was  developed.  The  parameters  obtained  from  the  identification 
procedure  are,  in  general,  the  coefficients  of  the  numerator  and 
denominator  polynomials  of  the  transfer  function  transformed  to  a 
V-plane.  A transform  associated  with  each  operator  used  was  defined 
such  that  the  transfer  function  can  be  transformed  back  from  the 
V-plane  to  the  S-plane.  From  the  general  class  of  operators  used,  two 
subgroups  of  operators  were  emphasized.  One  group  is  represented  by 
all  the  operators  preserving  the  S-plane  coefficients.  The  other  is 
the  group  of  operators  defined  in  the  discrete  z-domain. 

It  was  demonstrated  that  parametric  techniques  from  the  published 
literature,  are  special  cases  of  the  general  model  associated  with  a 
particulor  approximation  technique.  It  is  emphasized  in  this  disser- 
tation that  vhere  are  two  separate  stages  at  which  the  identification 
is  carried  out.  One  is  the  generation  of  the  basis  set  via  the 
chosen  linear  operator.  The  other  is  the  choice  of  an  approximation 
scheme. 

The  use  of  different  operators  and  different  approximation  schemes 
was  demonstrated  through  a variety  of  examples.  It  was  pointed  out 
that  some  operators  might  perform  tatter  than  others  in  a given  problem. 
In  general,  the  performance  of  an  operator  is  dependent  upon  the  data 
provided,  although  some  operators  are  shown  to  be  inherently  "bad 
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performers 

Having  in  mind  line  models  of  spectral  estimates,  there  is  a 
tendency  to  require  very  long  records  of  data  such  that  high  fre- 
quency "resolution1'  be  obtained.  This  is  true  for  line  spectra  and  is 
attributed  to  the  Sine  function  behavior  of  the  correlation  coefficient 
between  two  sinusoids  (i.e.,  the  correlation  function  decays  as  t ^). 
For  decaying  signals  like  exponentials,  this  conclusion  is  false.  It 
was  shown  that  for  exponentials  the  correlation  coefficient  asymptotes 
to  a value  defined  by  the  exponentials'  parameters.  Hence,  the  magni- 
tude of  the  correlation  coefficient  cannot  be  made  arbitrarily  small 
by  increasing  the  record  length.  A "long-record"  was  developed  based 
upon  the  decaying  signal  parameters.  This  "long-record"  is  of  finite 
length  and  any  increase  of  that  record  will  add  negligible  improve- 
ment to  the  identification.  Also,  on  the  basis  of  the  "long  record"  a 
rough  estimate  for  the  system's  order  was  developed. 

The  "growing  tree"  was  introduced  as  another  aid  in  determining 
the  system's  order.  The  potential  benefit  of  the  "growing  tree"  when 
applied  to  time  domain  synthesis  problems  was  also  pointed  out. 

Application  of  the  identification  model  to  slowly  time-variant 
systems  and  weakly  nonlinear  systems  was  discussed  and  demons :rated 
through  numerical  examples. 

5.2  Areas  for  Future  Research 

The  vast  bank  of  operators  should  be  investigated  and  classified 
on  the  basis  of  their  performance  when  applied  to  different  types  of 
identification  problems. 


The  applicability  of  the  identification  model  to  the  extended  class 


of  slowly  time-variant  systems  and  weakly  nonlinear  systems  should  be 
further  investigated. 


The  "growing  tree"  concept  appears  to  be  a promising  tool.  Its 
behavior  should  be  analyzed  on  firmer  mathematical  grounds.  From 
the  different  examples  worked  out,  it  seems  that  a signal’s  rate  of 
decay  and  its  bandwidth  are  strongly  related  to  the  region  in  the 
S-plane  occupied  by  the  tree.  This  relation  should  be  examined. 


Appendix  A 

AN  INTRODUCTION  TO  SOME  KNOWN  APPROXIMATION  TECHNIQUES 


-1.  Methods  Involving  a Known  Set  of  Basis  Functions 


Let  (g^Ct)},  i = 1,2,  ...,N,  denote  a given  set  of  basis  functions. 
Assume  that  the  real  time  function  x(t)  is  to  be  approximated  over  the 
closed  time  interval  [T^.T^]  by  the  finite  sum 


x(t)  * l Kg  (t). 
i-1 


(A- 1-1) 


The  problem  is  to  determine  the  N coefficients  K^,  i = 1,2,...,N,  so 
as  to  obtain  a reasonable  approximation.  Different  approaches  to  the 
approximation  will,  of  course,  result  in  different  sets  of  coefficients. 

Basically,  what  is  needed  is  a set  of  N linearly  independent 
equations  involving  the  N unknown  expansion  coefficients.  One  approach 
is  to  require  the  approximation  to  equal  x(t)  at  N different  instants 
of  time  t j , j = 1,2,...,N,  where  T^  t..  T2.  The  resulting  approxi- 

mation is  written  as 


*I(t)  - 1 UjgjCt) 

i=l 


(A-l-2) 


where  the  coefficients  are  determined  from  the  N linear  equations 


x^.(t^.)  * x(t^)  , J “ 1,2, .. . ,N. 


(A- 1-3) 


The  subscript  I in  (A-l-2)  is  used  to  indicate  that  the  approximation 
can  be  interpreted  as  an  interpolation  scheme  for  obtaining  approximate 
values  of  x(t)  at  instants  different  from  t ^ . This  approach  has  the 
serious  disadvantage  that  all  time  functions  having  the  same  sample 
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values  at  the  N instants  of  time  result  in  the  same  approximation.  There- 
fore, there  is  no  control  over  the  interpolation  error  at  the  nonsampling 
instants. 

The  approximation  may  be  improved  by  utilizing  M sample  points 
where  M > N.  We  can  no  longer  require  zero  error  at  the  sampling  in- 
stants. Nevertheless,  a "better"  overall  fit  can  be  obtained  according 
to  some  specified  criterion.  Using  a least-square  approach,  we  write 
the  approximation  as 

xT_(t)  = l L g (t)  (A- 1-4) 

Lb  i=l 

where  the  coefficients  L^,  i = 1,2,...,N  are  found  by  minimizing  the  sum 
of  the  squared  errors. 

A r ~ 2 

e = Z tx(t.)  - x (t  )]  . (A-l-5) 

IS  j 1 

j=l 

Carrying  out  the  minimization  with  respect  to  the  coefficients  results 
in  the  N linear  equations 


9e 


LS 


M 


N 


9L. 


-2  l [x(t  ) - l L g (t  )]g  (t  ) = 0,  k - 1,2 N. 

j=l  3 i=l  3 3 

(A-l-6) 


Rearranging  terms,  we  obtain 

M N M 

I I 8k^tj^8i^tj^Li  = ^ k ” 1,2,.. .,N.  (A-l-7) 

j 1 i 1 11 


When  M » N,  (A-l-7)  reduces  to  (A-l-3). 


A-3 


The  sampling  ratio  is  defined  to  be 


„ M 
R - — . 
N 


(A- 1-8) 


When  this  ratio  is  greater  than  unity,  the  problem  is  said  to  be  over- 
deternined.  Clearly,  the  average  distance  between  the  sampling  instants 
is  reduced  by  a factor  of  R in  the  overdetermined  case.  In  this  sense, 
better  control  of  the  interpolated  values  is  obtained  due  to  the 
higher  density  of  sampling  points.  However,  as  before,  all  time  func- 
tions having  the  same  sample  values  at  the  M time  instant  result  in  the 
same  approximation.  Ideally,  the  larger  the  value  of  R,  the  better. 

On  the  other  hand,  the  computation  becomes  excessive  when  R is  made 
too  large.  Also,  the  equations  become  linearly  dependent  due  to 
numerical  roundoff  as  the  sampling  instants  become  too  close. 

Instead  of  restricting  attention  to  discrete  instants  of  time, 

we  now  generate  approximations  by  considering  the  continuous  interval 

[TrT2].  The  integrated-squared-error  is  defined  to  be 
T„ 


ISE 


[x  (t ) - dt 


(A- 1-9) 


where 


xISE^t^  £ A^g^(t) 


(A- 1-10) 


The  coefficients  A , i = 1,2,...,N,  are  obtained  by  minimizing  e . 

1 lOL 

The  minimization  results  in 

T T 

2 2 

N r 

l t 8k(t)g1(t)dtJAi  = j gk(t)x(t)dt,  k-1,2,. . , ,N  . (A-l 

t -r 


"ID 


In  general,  it  is  not  possible  to  solve  for  each  coefficient  indi- 
vidually because  (A-l-11)  represents  N simultaneous  linear  equations. 
The  situation  is  greatly  simplified  when  the  basis  functions  are 
orthonormal. 

For  emphasis,  we  denote  a set  of  basis  functions  which  are  ortho- 
normal over  the  interval  [T^T^  by  {iJ/J  and  write  tb;  approximation  as 


;OISE(t)  ' Vl(t)- 


(A- 1-12) 


To  simplify  matters,  the  basis  functions  are  assumed  to  be  real.  The 


orthonormal  relationship  is  then  written  as 


- j ^i(t)^(t)dt  = 6.j 


(A-l-13) 


where 


1 i = j 

o i i j 


(A-l-14) 


Minimizing  the  integrated-squared-error  for  the  orthonormal  basis 


representation,  the  equation  for  the  coefficients  a.  become 


z 

ai  j *i 


(t)x(t)dt,  i = 1,2,... ,N 


(A-l-15) 


The  solutions  for  the  coefficients  are  seen  to  be  uncoupled  for  this 
special  case. 

As  with  the  previous  methods,  several  different  time  functions 


can  result  in  the  same  approximation.  For  example,  let  n(t)  be 
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orthogonal  to  {g^}  over  the  interval  [T^,T?].  It  then  follows  that 

the  integrated-squared-error  approximations  of  x(t)  and  y(t)  = x(t)  + n(t) 

are  identical.  However,  the  integrated-squared-errors  e (x)  and 

I SE 

elSE(y)  be  different  for  the  two  waveforms. 

Another  approach  to  the  approximation  problem  involves  considera- 
tion of  the  derivatives  and/or  integrals  of  x(t).  Let  the  n*’*'  derivative 
of  x(t)  be  denoted  by 


x(n)(t)  = d%(0  (A-l-16) 

dtn 

where  n is  a nonnegative  integer.  The  first  N terms  of  a Taylor  series 
expansion  about  the  point  t = T results  in  the  approximation 

x (t)  = [ x(£)(T)  T)  (A-l-17) 

P f=0 


where  the  subscript  p indicates  that  the  basis  functions  consist  of 
polynomials  in  t , the  error  is  given  by  the  remainder 


R^Ct)  = x ( t ) - xp(t) 


1 

(N-l)! 


t 


(t-T) 


N-l  (N) 
x 


J 

T 


(T)dT 


(A-l-18) 


Proof  of  this  follows  directly  by  evaluating  the  integral  in  (A-l-18) 
through  N successive  integrations  by  parts. 

To  obtain  an  accurate  approximation  of  x(t)  in  the  vicinity  about 
t = T when  using  a general  set  of  basis  functions  {g^},  we  require  that 
the  approximation  and  its  first  N-l  derivatives  equal  the  function  and 
its  first  N-l  derivatives,  respectively,  at  t = T.  Denote  the  approxi- 
mation by 


ll 


(A-l-19) 
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XED(C)  = l Ei8l(t) 
i»l 


where  the  subscript  ED  represents  equality  of  derivatives  at  time  T. 
Hence,  the  N coeffK  ients  E . , i = 1,2,...,N,  are  determined  from  the 
following  set  of  N linear  equations: 


>'  ' ' (T)  . I - 0,1 N-l 


(A-l-20) 


It  fol 


m.v.  be  expr  essed  as 

x<*>  (T) 

ED  U ’ l\ 


(A- 1-21) 


w t h*-  first  N terms  in  the  Taylor  series  expan- 
sion for  « (t).  B analogy  with  (A-l-18), 


ED 


(t)  " Xp(t)  = ~(N-~1)T  j(t;-T)N  1Xgo)(T)dT  . (A-l-22) 


Subtraction  of  (A-l-22)  from  (A-l-18)  yields  the  instantaneous  error 
expression 


£E(t)  x(t)  ~ XED(t)  = (N-l)! 


(A- 1-23) 


Let  the  interval  [T  ,T  ] be  contained  within  the  interval  [T  ,T  1. 

3 4 12 

An  approximation  to  x(t)  can  also  be  obtained  by  requiring  that  the 
approximation  and  its  first  (N-l)  integrals  equal  the  function  and 
its  first  (N-l)  integrals,  respectively,  over  [T^»T^].  The  n^-order 
integral  of  x(t)  i<  given  by 


t X3  X2 

x^  n^(t)  = ...  x(X  )dX  dX  ...  dX 

J 112  7t 

*3  3 3. 


(A-l-24) 


(t  - T ) 1 
(n-1) ! 


x(x)dT. 


The  approximation  is  denoted  by 


XEI(t)  = ^ Figi(t) 

t"L  i=l  1 1 


(A-l-25) 


where  the  subscript  El  indicates  equality  of  integrals  over  [T,,T.]. 

J 4 

The  coefficients  F^,  i = 1,2,...,N,  are  determined  from  the  N linear 
equat ions 


x{  £)(T  ) = x(  £)(T,),  i = 0,1,..., N-1. 
El  4 4 


(A-l-26) 


The  two  previous  approximations  may  be  combined  to  yield  yet 
another  approximation.  This  is  accomplished  by  requiring  the  first  K 
integrals  from  (A-l-26)  and  the  first  N-K  derivatives  from  (A-l-20)  to 
be  equal.  The  approximation  is  now  denoted  by 


Wt}  = Jx  Hi8i(t) 


(A-l-27) 


where  the  subscript  EID  indicates  equality  of  both  integrals  and  deriva- 
tives. For  simplicity  we  set  T^  = T.  The  coefficients  , i = 1,2,...,N, 
are  then  determined  from  the  N linear  equations 

x££J(T)  = x (£)  (T)  , Z « -K -1,0,1 N-K-l  . (A-l-28) 
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The  approximation  in  (A-l-27)  can  be  further  generalized  by  over- 


determining the  problem.  Assume  the  number  of  sampling  points  is 
given  by  M and  that  a total  of  L derivatives  and/or  integrals  of  x(t) 
are  to  be  considered.  Let  the  approximation  be  expressed  as 


XGLS 


(t) 


N 

l Qiei(t) 

i=l 


(A-J-29) 


where  the  subscript  GLS  indicates  that  the  coefficients  are  to  be 
determined  from  a generalized  least-square  approach.  Taking  into  con- 
sideration the  first  K integrals  and  the  first  L-K  derivatives  of 
x(t)  (including  the  zero-order  derivative),  the  sum  of  the  squared 
errors  becomes 


L-K-l  M f 

' Jk  ‘V 


.,2 

XGLS  ^ j ’ 


(A- 1-30) 


where  the  product  of  L and  M is  chosen  to  be  greater  than  or  equal  to  N. 
Notice  that  (A-l-30)  reduces  to  (A-l-5)  when  L = 1 and  K = 0.  Also, 
for  L - N and  M = 1,  the  problem  is  no  longer  overdetermined  and  determi- 
nation of  the  coefficients  by  the  generalized  least-squares  approach 


reduces  to  (A-l-28)  with  the  result  that  the  error  _ is  identically 

bLb 

zero.  Minimization  of  the  error  in  (A-l-30)  with  respect  to  the  coeffi- 
cients Q^,  1 = 1,2,...,N,  results  in  the  N linear  equations 


L-K-l  M N 

l l l 

A— K j-1  i=l 


gW(t  )g(X,)(t  )Q 
sk  v j/Bi  v j'^i 


L-K-l  M 

l I 

A— K j-1 


k = 1,2,. ...N 
(A-l-31) 


For  L = 1,  K = 0 and  M = N,  the  problem  becomes  the  interpolation  prob- 
lem considered  in  (A-l-2)  and  (A-l-31)  reduces  to  (A-l-3). 


A 


Similarly,  the  integrated  squared  error  approach  can  be  generalized 

by  involving  the  first  K integrals  and  L-K  derivatives  of  x(t).  The 

approximation  is  now  denoted  by 

N 

xGisE(t)  = l B.g.(t)  (A- 1-32) 


where  the  subscript  GISE  stands  for  generalized  integrated-squared-error. 
The  generalized  integrated-squared-error  is  defined  io  be 


eGISE 


L-K-l 


1 

£=-K 


x(£)(t)]2dt 


(A-l-33) 


Performing  the  minimization  of  (A-l-33)  with  respect  to  the  coefficients 
EL,  i = 1,2,...,N,  yields  the  N linear  equations 


L-K-l  N 

l l [ 

£— K i=l 


2 (£)  (£)  '2 

gi't;(t)gJJt;(t)dt]B.  = l 

k 1 £— K 


g^£)(t)x(£)(t)dt,  k = 1,2,...  , 
(A- 1-34) 


The  derivatives  and/or  integrals  of  x(t)  are  no  longer  involved  when 
L = 1,  K = 0 and  (A- 1-34)  becomes  (A-l-11). 

As  a final  remark,  we  point  out  that  the  sampling  instants 
appearing  in  the  set  of  N equations  of  (A-l-20)  can  be  allowed  to  vary 
from  one  equation  to  another.  The  same  statement  applies  to  the  sets 
of  N linear  equations  appearing  in  (A-l-26),  (A-l-28),  and  (A-l-31). 
This  allows  for  greater  flexibility  in  generating  the  approximations. 

At  this  point  we  illustrate  the  previous  discussion  with  a simple 
example.  The  objective  is  to  demonstrate  the  varying  amounts  of  compu- 
tational effort  required  and  the  different  waveforms  which  result  using 
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the  various  methods  of  approximation.  Following  the  example,  we  discuss 
how  the  different  approximations  are  related  and  show  that  a natural 
criteria  for  their  comparison  is  the  integrated  square  error. 

Example  A- 1-1 

Given  the  finite  set  of  basis  functions 

{gt}  = {l,t,t2}  (A-l-35) 

approximate  the  function 

x(t)  = cos  2irt  (A-l-36) 

over  the  interval  = [0,1]. 

Case  1 (interpolation  scheme) 

The  approximation  is  given  by  (A-l-2)  where  the  coefficients  are 
obtained  from  (A- 1-3).  The  sampling  instants,  which  are  equal  in 
number  to  the  number  of  coefficients  to  be  determ. ned,  are  arbitrarily 
cbsen  to  be  t^  = 0,  1/2,  1.  Sample  values  of  both  x(t)  and  the  basis 
functions  are  given  in  Table  A- 1-1.  In  matrix  form,  (A-l-3)  becomes 


Table  A- 1-1 


/~s 

V-/ 

X 

gl 

g2 

g3 

0 

1 

1 

0 

0 

1/2 

-1 

1 

1/2 

1/4 

1 

1 



1 



1 

1 

1 

1 

1 


0 

1/2 

1 


0 

1/4 

1 


(A- 1-37) 


Solution  of  (A-l-37)  for  the  coefficients  yields  the  approximation 


~ 2 

XjCt)  = 1 - 8t  + 8t  (A- 1-38) 

x(t)  and  x (t)  are  shown  in  Fig.  A- 1-1. 

Ic 


o 


/\ 

Fig.  A-l-1.  The  function  x(t)  and  its  approximation  Xj(t). 

Case  2 (least-squares  approach) 

The  approximation  is  given  by  (A-l-4)  where  the  coefficients  are 
determined  from  (A- 1-7).  We  arbitrarily  choose  a sampling  ratio  of 
3 so  that  M = 9 and  N = 3.  The  9 sampling  instants  and  the  corresponding 
sample  values  of  both  x(t)  and  the  basis  functions  are  shown  in 
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"able  A-l-2.  Let  the  9x3  matrix  formed  from  the  columns  in  Table 

A-l-2  be  denoted  by  [G].  In  terms  of  matrices,  (A-l-7)may  be  expressed 


as 


T 

[G]  [G] 

L1 

L2 

II 

o 

' H 

xCt^ 

L3 

X(tg) 

(A-l-39) 


Table  A-l-2 


t . 

J 

x(t.) 

J 

81 

g2 

g3 

0 

1 

1 

0 

0 

1/6 

1/2 

1 

1/6 

1/36 

1/4 

0 

1 

1/4 

1/16 

1/3 

-1/2 

1 

1/3 

1/9 

1/2 

-1 

1 

1/2 

1/4 

2/3 

-1/2 

1 

2/3 

4/9 

3/4 

0 

1 

3/4 

9/16 

5/6 

1/2 

1 

5/6 

25/36 

1 

1 

1 

1 

1 

Solution  of  (A-l-39)  for  the  coefficients  results  in  the  approximation 

xLS(t)  = 1.1668  - 7 . 0524t  + 7.0524t2.  (A-l-40) 

A 

x(t)  and  x^<,(t)  are  shown  in  Fig.  A-l-2. 

Case  3 (integrated  squared  error  approach) 

The  set  of  basis  functions  (g.),  given  by  (A-l-35),  does  not  con- 
stitute an  orthonormal  set.  To  simplify  evaluation  of  the  coefficients, 
the  Gram-Schmidt  orthonormalization  procedure  was  applied  to  {g^} 
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x(t)  and  x0ISE(t)  are  shown  in  Fig.  A-l-3. 

t 


o 

o 


(X 1 0-1 J 


Fig.  A-l-3.  The  function  x(t)  and  its  approximation  x (t). 

OISE 

Case  4 (equality  of  derivatives  approach) 

The  approximation  is  given  by  (A-l-19)  where  the  coefficients  are 
obtained  from  (A- 1-20).  As  is  apparent  from  (A-l-20),  the  basis  func- 
tions and  their  derivatives  are  needed  for  evaluation  of  the  coeffi- 
cients. Since  three  basis  functions  are  included  in  {g^},  N=3. 
Therefore,  (A-l-20)  consists  of  three  simultaneous  equations  for  which 


l = 0,1,2.  The  corresponding  basis  functions  are  given  in  Table  A-l-3. 
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Table  A-l-3 


£ 

g(<0 

g(*} 

s2 

gW 

0 

i 

t 

2 

t 

1 

0 

1 

2t 

2 

0 

0 

2 

In  addition,  the  signal  to  be  approximated  and  its  first  two  derivatives 
are 

x(0) (t)  = cos  2nt 

x^(t)  = — 2tt  sin  2nt  (A-l-43) 

(2)  2 

xv  y(t)  = -4tt  cos  2irt. 

Applying  (A-l-20)  at  the  time  instant  t = T » 1/2,  solving  for  the  co- 
efficients E^,  and  substituting  into  (A-l-19)  results  in 

2 

/n  w - 2 2 2 2 

XED(t)  = " 2 " 2lT  C + 2lT  • (A- 1-44) 

A A 

x(t)  and  xED(t)  are  shown  in  Fig.  A-l-4.  Note  the  close  fit  of  xED(t) 
to  x(t)  around  t = 1/2. 

Case  5 (equality  of  integrals  approach) 

The  approximation  is  now  given  by  (A-l-25)  where  the  coeffl  ients 
are  determined  from  (A-l-26).  As  in  Case  4,  three  simultaneois  equa- 
tions are  needed  since  N = 3.  Integrating  over  the  interval  [1,'  ],  the 
basis  functions  corresponding  to  £ = 0,1,2  are  shown  in  Table  A- 1-4. 


Solving  (A-l-26)  with  T “ 1 and  T^  = 0,  the  approximation  becomes 


xEI(t)  = 1 - 6t  + 6t  . 


x(t)  and  x£I(t)  are  shown  in  Fig.  A-l-5. 


(A-l-46) 


0.00  2.50 

(XI 0-1) 


5.00 


7.50 


10.  00 


Fig.  A-l-5.  The  function  x(t)  and  its  approximation  x (t). 
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Case  6 (equality  of  integrals  and  derivatives  approach) 

For  this  case  the  approximation  is  given  by  (A-l-27)  where  the  co- 


efficients are  evaluated  from  (A- 1-28).  Once  again,  N = 3.  To  obtain 
these  three  simultaneous  equations,  we  let  K - 1 in  (A-l-28).  Choosing 
T = 1/2  and  evaluating  the  integral  over  the  interval  (1,1/2),  the 


basis  functions  evaluated  at  T = 1/2  are  given  in  Table  A-l-5. 


Table  A-l-5 


i 

gj S) (1/2) 

g2(£)(l/2) 

g^}(l/2) 

-1 

-1/2 

-3/8 

-7/24 

0 

1 

1/2 

1/4 

1 1 0 

1 

1 

Similarly , 

k('1)(1/2)  = 0,  x(0)  (1/2)  = -1,  x(1)(l/2)  = 0 . (A- 1-47) 

This  results  in  the  approximation 

xEID(t)  = 2 “ 12t  + 12t2.  (A-l-48) 

/\ 

x(t)  and  xEID(t)  are  shown  in  Fig.  A-l-6. 


(X  1 0-1) 

Fig.  A-l-6.  The  function  x(t)  and  its  approximation  x (t). 

E1D 


A-19 


Case  7 (gener alized  least-square  approach) 

The  approximation  to  x(t)  is  given  by  (A-l-29)  where  the  coefficients 
are  obtained  from  (A-l-31).  We  present  two  different  illustrations. 

In  the  first  a single  sampling  point  (M  = 1)  is  utilized  whereas  9 dif- 
ferent sampling  points  (M=  9)  are  used  in  the  second.  In  both  illus- 
trations, we  consider  the  function  and  4 of  its  integrals  over  the 
interval  (l,t). 

Part  A:  In  this  illustration,  N = 3,  M = 1,  K - 4,  and  L = 5.  The  basis 

(£) 

functions  g^  for  £ = 0,-1, -2  and  i * 1,2,3  are  given  in  Table  A-l-4. 

The  functions  corresponding  to  £ = -3,-4  are  listed  in  Table  A-l-6.  The 
function  being  approximated  and  its  first  two  integrals  are  given  in 


Table  A-l-6 


£ 

8f> 

8a) 

g2 

gW 

83 

-3 

(t-D3 

(t-1) 3(t+3) 

(t-l)3(t2+3t+6) 

6 

24 

60 

(t-1)4 

(t-l)4(t+4) 

(t-l)4(t2+4t+10) 

24 

120 

360 

(A-l-45) . In  addition,  xv  3\t)  and  x^  ^(t)  are 


X(~3)(t)  - (t-1)  - Sin  2TTt 

4tt  8tt 

x^  ^(t)  = (t-1)2 (1  - cos  2nt)  . 

8tt  16ir 

(£) 

The  sampling  instant  and  the  corresponding  sample  values  of  xv 

(£) 

and  g^  (t)  are  shown  in  Table  A- 1-7. 


(A- 1-49) 


(t) 


For  a specific  value  of  £,  let  the  1x3  matrices  formed  from  the  columns 

(£)  (£) 

be  denoted  by  [G  ].  In  terms  of  matrices,  (A- 1-31)  may  be  ex- 
pressed as 


{ l |G(I)]T[G<,!)])  Q,  - I [GW]V‘,<t) 


'l  0 


a>,T  aj, 


(A- 1-50) 


Solution  of  (A-l-50)  for  the  coefficients  results  in  the  approximation 


x (t)  = 1.00013  - 6. 12286t  + 6.18055t  . 

oLo 


(A- 1-51) 


x(t)  and  x (t)  are  shown  in  Fig.  A-l-7. 
GLb 


0.00  2.50 

(X10~U 
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10.00 


Fig.  A-l-7.  The  function  x(t)  and  its  approximation  x„TC(t),(  =1) , 
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Part  B:  Once  again,  N = 3,  K = 4,  and  L = 5.  However,  now  rl = 9 since 
we  sample  at  9 different  time  instants.  The  sampling  instants  and  the 
corresponding  sample  values  of  x^  ^ (t)  and  gf^(t)  are  shown  in  Table 
A-l-8.  For  a specific  value  of  l , let  the  9x3  matrices  formed  from  the 
columns  g^  be  denoted  by  [Gv  ].  In  terms  of  matrices,  (A-l-31) 
becomes 


0 

{ l tG 

£.“-4 


<*>]T[gW 


■Q1 ' 

o 

Q2 

- 1 [C(t)|T 

[<*!j 

II 

1 

(A-l-52) 


Solution  of  (A-l-52)  for  the  coefficients  results  in  the  approximation 

XcLgCt)  = 1.17  - 7. 1014t  + 7.l089t2  . (A-l-53) 

A 

x(t)  and  x,.TC(t)  are  shown  in  Fig,  A-l-8. 


Fig.  A-l-8.  The  function  x(t)  and  its  approximation  x (t),  (M  = 9). 


A-  U 


Case  8 (generalized  integr ated-squ ared-error  approach) 

The  approximation  to  x(t)  ig  given  by  (A-l-32)  where  the  coeffi- 
cients IK  are  obtained  from  (A-l-34) . Since  the  coefficients  in  the 
(USE  approach  are  determined  by  minimizing  the  sum  of  errors  indicated 
by  (A- 1-33),  it  is  of  interest  to  investigate  the  magnitude  of  the 

individual  errors  in  the  sum.  The  values  of  the  coefficients  B are 

1 

most  influenced  by  the  need  to  minimize  the  larger  errors. 

For  this  case  the  approximation  is  generated  by  considering  the 
original  function,  its  derivative,  and  its  integral  (i.e.,  L = 3,  K = 1). 
To  assess  the  size  of  the  errors  likely  to  be  encountered  in  (A- 1-33) , 
we  first  carry  out  ISE  approximations  to  (1)  ^(t)  in  terms  of 

{g|  ^(t)},  (2)  x(t)  in  terms  of  (g^(t)},  and  (3)  x^(t)  in  terms  of 
(gj^  (t)}.  As  pointed  out  earlier,  the  (USE  approach  can  be  specialized 
to  the  ISE  approach  with  the  proper  choice  of  L and  K in  (A-l-34) . 


(1)  approximation  of  x^  ''ft)  in  terms  of  {g  j ^(t)}  (L=l,  K=l) 

The  inner  product  for  this  case  is  defined  to  be 

1 

,>a  f 


<y,  z: 


! y(t)z(t)dt  . 
' 0 


(A- 1-54) 


Hence,  the  inner  produi ts  required  by  (A-l-34)  are: 


Kg 

L 3 


.S, 


(-1) 

(-1), 

„ (-D 

(-1)  > 

11 

5 

r 

,R2 

<8l 

,g3  | 

|3 

24 

(-D 

2 

5 

= 24 

2 

15 

(-1) 

3 

,g2(-1)> 

l_2 

20 

7 

72 

(-1)  . (-1). 

k 


[<8i)  ’ ,x 


1 

, (-1)  (-l)s 

= (-  -*y> 
4tf2 

1 

1 

<g2  *x 

2 

4tt2-6 

12it2 

(U —1—55) 


0|U> 


In  terms  of  the  matrices  of  (A-l-55),  (A-l-34)  becomes 


.(-1)  „(-D 


(-1)  (-D. 


[<g£  *',g } X'>]  B2  = [<g^  A',xv  x'>]  . 


(A-l-56) 


The  approximation  is  given  by 


xISE)(t)  = 1-787(t"1)  " 10.332[^(t2-D]  + 10. 211  (t3-l)] 


= -.02427  + 1. 787t  - 5.166t2  + 3.404t3 


(A-l-57) 


(-1)  ~(-l) 

x and  Xjggit)  are  shown  in  Fig.  A-l-9. 


'o.oo  2.50  5.00  7.50  10.00 

txio-n 


(-13  '-(-13 

Fig.  A-l-9.  The  function  x (t3  and  its  approximation  x.  (t). 
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2)  approximation  of  x(t)  in  terms  of  {g^(t)}  (L*l,  K=0) 

The  inner  products  required  by  (A-l-34)  are  now  given  by: 


’<81 

<gi,g2> 

<g1.g3> " 

1 

1 

2 

[<gk,gi>]  = 

<g2 ,81> 

<g2»g2> 

<g2«g3> 

= 

1 

2 

1 

3 

<g3,gl> 

<'g3»g2> 

*^83  *g3> 

1 

3 

1 

4 

<g1>x> 

’ 0 

[<gJc»x>]  = 

<g2>x> 

<g3»x> 

= (72} 

2tt 

0 

1 

(A-l-57) 


The  matrix  equation  for  the  coefficients  B.  is  similar  in  form  to 
(A-l-56)  and  results  in  the  approximation 


XISE 


90  90 

Tr2  "h-2 


2 

t 


= 1.51982  - 9. 11891t  + 9.11891t2. 


(A-l-58) 


This  is  identical  to  the  approximation  that  resulted  in  Case  3.  (See 
Fig.  A-l-3) . 


3)  approximation  of  x^^(t)  in  terms  of 


(L=l ,K=-1) 


The  inner  products  required  by  (A-l-34)  are  now  given  by 


- 


[<g^,gj1)>  <£>,£>>  <g$1),41)> 


Kg^.x^J 


U) 

U) 

’1 

,gl 

(1) 

B(1) 

7 

,81 

(1/ 

(1) 

■3 

,81 

(1) 

(1) 

1 

.x 

,(1) 

(1) 

*2 

.X 

(1) 

(1) 

•3 

.X 

for 

the 

^ .g2  > <82  »g3 


"~g2  »82  g2  .83 


0 0 0 


Oil 


0 1 5 


(A-l-59) 


and  results  in  the  approximation 


x^'(t)  = c(0)  - 6(1)  + 6(2t) 


(A-l-60) 


-6  + 12t  . 


(1)  ''(l) 

where  c is  an  arbitrary  constant.  x (t)  and  X^gE(t)  are  shown  in 


Fig.  A-l-10. 


ISE 
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Fig.  A-l-10.  The  function  x^(t)  and  its  approximation  x5i^.(t) 
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Having  obtained  these  three  approximations,  note  that  the  error  becomes 

y'(-l)  ''(l) 

progressively  larger  as  one  proceeds  from  x'  (t)  to  x (t). 

ISE  IS  E 

We  now  approximate  x(t)  by  means  of  the  GISE  approach.  However, 
to  emphasize  the  importance  of  the  various  terms  in  (A-l-33),  we  first 
consider  the  situation  for  which  L = 2,  K = 1.  With  these  values  of 
L and  K only  x(t)  and  x^  "^(t)  are  involved  in  the  approximation.  We 
then  repeat  ‘ he  procedure  for  L = 3,  K = 1.  Then  x(t),  x^  ^ (t)  , and 


xv  (t)  are  all  involved. 

Equation  (A-l-34)  may  be  expressed  in  matrix  form  as 

T <),gf)>]  5 =T  [<8f),xw>] 

ft— K k 1 „ ft— K 

B3 


(A-l-61) 


For  L * K=  1,  the  approximation  is 

X-- _ ( t ) = 1.522113  - 9. 132479t  + 9.132466t' 
GISE 


(A- 1-62) 


Note  that  the  coefficients  in  (A-l-62)  fall  between  those  in  (A-l-57) 

/\ 

and  (A-l-58) . However,  since  the  error  associated  with  xTC,„(t)  is 

IbE 

-(-1) 

greater  than  that  associated  with  x (t),  the  coefficients  fall  closer 

IbE 

to  those  in  (A-l-58).  For  L = 3,  K = 1,  the  approximation  is 


xC[SF(t)  » 1.022  - 6 . 0523t  + 6.0525t 


(A-l-63) 


The  coefficients  in  (A-l-63)  now  fall  between  those  in  (A-l-57)  and 

(A-l-60).  However,  since  the  error  associated  with  *jgg(t)  is  much 

larger  than  either  of  the  errors  associated  with  x^~^(t)  and  xTC„(t), 

IbE  ISE 

the  coefficients  fall  very  close  to  those  in  (A-l-60).  Plots  of  x(t), 
XGISE(t)|L-2,K=l*  and  XGISE(t)lL=3,K-l  are  shoWn  in  Flg*  A"1"11- 
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Fig.  A- 1-11.  The  function  x(t)  and  its  approxinations 


XGISE(t) k=2,K=l  ’ XGISE^t) Il=3,K=1' 

The  relationships  between  the  various  approximations  are  indicated 
by  the  block  diagram  presented  in  Fig.  A-l-12.  The  equation  used 
to  solve  for  the  coefficients  is  refer. ed  to  in  each  block.  The  directed 
paths  between  the  blocks  serve  to  indicat.  how  one  approach  may  be 
transformed  into  another. 

Note  that  the  block  diagram  in  Fig.  A-l-12  subdivides  the  various 
approximations  into  those  that  process  continuous  waveforms  versus 


numerical  data.  In  general,  as  the  number  of  sampling  instants  is 
increased  (l.e.,  R , the  solution  for  the  discrete  case  converges 

to  that  of  the  continuous  case.  This  is  illustrated  in  Fig.  A-l-13 
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where  L^(R)  denotes  the  i1,1  least-squares  coefficient  of  (A-l-7)  for 
a specified  value  of  R and  A^  denotes  the  i^f'  integrated-squared-error 
coefficient  of  (A-l-11).  In  our  example,  as  R was  increased  with  a 
fixed  value  of  N = 3,  the  smallest  magnitudes  of  the  coefficients  L^, 
L^,  and  all  occurred  for  M = 4.  As  a result,  L.(R)  was  normalized 
by  (4/3)  in  order  to  compress  the  plot.  The  LSE  coefficients  nicely 
asymptote  to  the  ISE  coefficients.  The  convergence  is  even  more 


dramatically  illustrated  by  examining  the  errors.  Because  the  wave- 
form being  approximated  is  known  in  our  example,  it  is  possible  to 
compute  the  integrated-squared  error  for  both  the  discrete  and  continuous 
cases.  These  are  shown  in  Fig.  A-l-14.  The  asymptote  is  determined 
from  Ejgg>  the  integrated-squared  error  for  the  continuous  case.  The 
integrated-squared  error  for  the  discrete  case,  for  a specified  value 
of  R,  is  denoted  by  For  convenience,  all  errors  are  normalized 

by  that  occurring  for  the  interpolation  case  (i.e.,  M = N = 3).  Note 
the  abrupt  convergence.  The  plot  reveals  that  a value  of  R equal  to  10 
is  more  than  adequate  for  the  discrete  case  to  approximate  the  continuous 


A summary  of  the  waveforms  resulting  from  the  various  approximations 
is  given  in  Fig.  A-l-15.  The  associated  integrated-squared  error  is  in- 
dicated on  each  plot.  By  examination  of  Case  4,  it  is  obvious  that  an 

"interval  type"  of  error  criterion,  such  as  e , is  inappropriate  for 

ISE 


"Taylor  type"  of  approximation. 


.20 
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R“M/N 

Fig.  A-l-13.  Normalized  approximation  coefficients  Vs.  sampling 
ratio. 
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Fig.  A-l-14.  Normalized  integrated  squared  error  Vs.  sampling  ratio. 
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A.  2 Convergence  as  R •»  00  with  a known  exponential  basis. 

As  a second  example,  the  rate  of  convergence  of  the  coefficients 
and  error  in  the  discrete  problem  to  those  of  the  continuous  problem 
is  investigated  when  the  chosen  basis  is  of  an  exponential  type.  For 
simplicity,  the  coefficients  are  restricted  to  be  real.  This  con- 
straint is  met  by  combining  complex  exponentials  with  their  conjugates 
in  order  to  obtain  decaying  sinusoids. 


Example  A-2-1 

Given  the  finite  set  of  basis  functions 

r , r - (l-jTT)t  ~ ( 1+ j TT  ) t -(2-j3*)t  - ( 2+ j 37T  ) t -(3-j  57T  ) t -(3+j5TT)t 
1 8± 1 = le  .e  , e ,e  ,e  ,e 


where,  j * /-l,  approximate  the  function 


x(t)  = 


1 ; 0 < t < 1 


0 , otherwise 


(A-2-1) 


(A-2-2) 


over  the  interval  [T^,  T^]  = [0,2], 

Since  the  coefficients  in  (A-l-7)  are  constrained  to  have  a 
zero  imaginary  part  and  since  x(t)  is  real,  the  complex  conjugate  basis 
functions  represented  in  (A-2-1)  can  be  combined  into  a new  set  of  real 
basis  functions  composed  of  decaying  sinusoids.  Consequently  an  equi- 
valent set  of  basis  functions  is  given  by 

2 1 -2t  — 3 1 

{g  }=  (e  cos  irt,  e sin  TTt,  e cos  3iTt,e  sin  3irt,e  cos  5iTt} 
Si  (A-2-3) 

where  e ^tsin  5irt  has  been  intentionally  omitted  in  order  to  increase 
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the  approximation  error. 

As  in  example  A- 1-1,  the  waveform  being  approximated  is  known. 
Hence,  computation  of  eISg(R)  is  possible  for  the  discrete  case.  The 
normalized  least-square  coefficients  L^CR),  i = 1,2, 3, 4, 5 are  plotted 
in  Figures  A-2-1  through  A-2-5,  along  with  the  asymptotes  which  are 
given  by  the  integrated-squared-error  coefficients  A^ . The  integrated- 
squared  error  for  the  discrete  case,  as  a function  of  the  sampling 
ratio,  is  shown  in  Fig.  A-2-6  where  the  asymptote  i-s  the  ISE  of  the 
continuous  case.  As  in  example  A- 1-1,  notice  that  R = 10  is  more  than 
adequate  for  the  discrete  case  to  approximate  the  continuous  case. 

Plots  of  x(t)  and  its  integrated-squared-error  approximation  are  shown 
in  Fig.  A-2-7. 

To  provide  additional  insight  into  the  convergence  of  the  least- 
squares  approximation  of  x(t),  as  given  in  (A-2-2) , the  approximation 
is  repeated  five  times  as  the  number  of  basis  functions  is  gradually 
increased  from  one  to  five.  In  particular,  the  following  sets  of  basis 
functions  were  used: 


(a)  {gsl} 

(b)  {gsl’gs2} 

(C)  {gsl’gs2’gs3} 


(A-2-4) 


Fig.  A-2-5.  Normalized  L-  coefficient  Vs.  sampling  ratio 


Fig.  A-2-6.  Normalized  integrated  squared  error  Vs.  sampling  ratio 
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f § 


Fig.  A-2-7.  x(t)  and  its  MISE  approximation. 

between  3 and  5 and  is  insensitive  to  the  number  of  basis  functions  used 
(i.e.,  insensitive  to  the  quality  of  the  approximation). 


A. 3 General  Discussion  of  Numerical  Examples 

From  the  theory  of  band-limited  signals,  it  is  well-known  that  one 
has  to  sample  a given  signal  at  least  at  the  Nyquist  rate  in  order  to 
achieve  a meaningful  reconstruction.  If  a signal  la  sampled  below  this 
rate,  an  aliasing  error  occurs.  This  type  of  error  can  be  interpreted 
in  terms  of  Riemana surfaces.  It  is  clear  that  a phase  variation  of 
larger  than  it  radians  between  consecutive  time  samples  for  any  of  the 
complex  exponential  waveform  components  will  locate  the  phase  of  that 
component  on  a Riemana surface  different  from  the  first.  If  no  informa- 
tion is  supplied  about  the  true  Riemanx surface  on  which  the  phase  is 


located,  then  an  error  will  occur  in  determination  of  that  component’s 
argument.  In  practise,  an  aliasing  error  always  occurs  since  physical 
signals  can  never  be  strictly  band  limited.  For  such  signals  the 
Nyquist  rate  is  usually  determined  by  specifying  a frequency  above  which 
the  energy  content  is  assumed  to  be  negligible.  The  aliasing  error 
is  then  a function  of  the  energy  ratio  between  the  in-band  and  out-of- 
band  components. 

Hence,  having  chosen  a suitable  set  of  basis  functions,  one  must 
assure  in  the  Interpolation  case  that  the  sampling  rate  is  chosen  to  be 
at  least  the  Nyquist  rate  if  a meaningful  reconstruction  is  to  be  ob- 
tained. In  the  preceding  examples,  it  was  shown  that  a sampling  ratio 
of  about  10  yields  discrete  approximations  having  an  approximation  error 
almost  as  small  as  that  of  the  continuous  MISE  approximation.  In  the 
Interpolation  problem  the  minimum  number  of  samples  required  equals 
the  number  of  known  basis  functions.  For  a band-limited  signal,  the 
Nyquist  rate  determines  the  minimum  number. of  samples  over  a given 
Interval.  This  implies  that,  in  general,  the  minimum  number  of  basis 
functions  for  a band-limited  signal  should  equal  the  Nyquist  rate  times 
the  length  of  the  sampling  interval  if  a meaningful  reconstruction  is 
to  be  obtained  in  the  interpolation  problem.  To  obtain  an  ISE(R)  error 
which  is  close  to  the  MISE  with  such  a set  of  basis  functions,  the 
waveform  should  be  oversampled  at  10  times  the  Nyquist  rate.  Hence, 

20  samples  per  cycle  of  the  highest  frequency  component  of  significance, 
speaking  energywise,  are  required.  If  the  ratio  of  highest  to  lowest 
frequency  is  100  and  the  approximation  is  carried  over  a full  cycle  of 
the  lowest  frequency,  2000  samples  are  needed.  Although  this  is  *• 
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extremely  large  number  of  samples,  one  should  remember  that  such  a 
ratio  of  frequencies  describes  a wide  band  signal.  Also,  the  value 
of  10  is  a conservative  choice  for  R.  For  such  wide-band  signals,  it 
may  be  desirable  to  pick  R as  low  as  3 to  5 (i.e.  600  to  1000  samples). 


Appendix  B [32] 


Exact  solution  to  the  synthesized  problem  of  Example  3-2-1. 
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SAMPLE 

TIME  (SEC) 

INSTANTANEOUS  FREQ  (Hz) 

t, 

.04 

141.47 

*, 

.08 

145.73 

t, 

.12 

152.37 

t. 

.16 

160.73 

t. 

20 

170.00 

t. 

2i 

179.27 

t, 

28 

187.63 

t. 

.32 

194.27 

t. 

.36 

198.53 

u 

.40 

200.00 

Ml 

.44 

198.53 

MI 

.48 

194.27 

*U 

.52 

187.63 

.56 

179.27 

** 

.60 

170.00 

* If 

.64 

160.73 
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